import sys
print(sys.version)
3.6.9 (default, Oct 8 2020, 12:12:24) [GCC 8.4.0]
from netCDF4 import Dataset, num2date
from datetime import datetime
import glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import matplotlib.dates as mdates
import calendar
import csv
from scipy.stats import describe
from scipy import interpolate
from tqdm import tqdm_notebook as tqdm
import oceansdb
from mpl_toolkits.basemap import Basemap
#import seawater as sw
#from seawater.library import T90conv
def get_first(d, varnames):
for var in varnames:
try:
return d.variables[var]
except Exception as e:
pass
raise KeyError(f"none of {varnames} in {d.variables.keys()}")
def load_netcdf(path = "Argo_South_60"):
# load_netcdf function input: path - default is "Argo_South_60"
files = glob.glob("data/" + path + "/**/*.nc", recursive=True)
# glob.glob - returns a list of all the files in a given folder with the extension 'nc'
# recursive=true - keeps on looking for all the files till the end, or the last file that exist
# creating lists working on objects
data = {
"lat": [],
"lon": [],
"datetime": [],
"pressure": [],
"temperature": [],
"salinity": []
}
datetimes = []
for f in tqdm(files):
# tqdm - make your loops show a progress meter
d = Dataset(f)
####filename = "{}".format(f)
# DATASET - Creating/Opening/Closing a netCDF file.
#This is the method used to open an existing netCDF file.
lat = get_first(d, ["LATITUDE", "lat"])[:]
mask = lat < -60
lon = get_first(d, ["LONGITUDE", "lon"])[:]
#1 The extend () method adds the specified list elements to the end of the current list.
#2 This method does not return any value but add the content to existing list.
juld = get_first(d, ["JULD", "time"])
units = juld.getncattr('units')
#1 getncattr - reads the attributes of given variable i.e units
dates = num2date(juld[:], units, "standard", only_use_cftime_datetimes=False)
#1 num2date - converts the format of the loaded datetime to standard format YYYY-MM-dd
try:
data["pressure"].extend(get_first(d, ["PRES_ADJUSTED", "Pressure", "depth"])[:][mask])
data["temperature"].extend(get_first(d, ["TEMP_ADJUSTED", "Temperature"])[:][mask])
data["lat"].extend(lat[mask])
data["lon"].extend(lon[mask])
try:
data["datetime"].extend(dates[mask])
except:
if len(np.flatnonzero(mask)):
data["datetime"].append(dates)
try:
data["salinity"].extend(get_first(d, ["PSAL_ADJUSTED", "Salinity"])[:][mask])
except:
data["salinity"].extend(np.full(len(mask), np.nan))
except Exception as e:
print(f"Skipping {f} due to {e}")
#1 array() - allows one to do operations on a list
#2 NumPy was imported as np
#3 If we want to do math on a homogeneous array of numeric data, then it is much better to use NumPy,
# which can automatically vectorize operations on complex multi-dimensional arrays
#1 array() - allows one to do operations on a list
#2 NumPy was imported as np
#3 If we want to do math on a homogeneous array of numeric data, then it is much better to use NumPy,
# which can automatically vectorize operations on complex multi-dimensional arrays
# k =
#df = pd.DataFrame(datetimes)
#print ("dataframe df[0]"); print(df[0])
#df = df.sort_values(by=[0]).reset_index()
#df = df.value_counts(by=[0]).reset_index()
#print("df sorted");print(df[0])
#return datetimes, df
#df.to_csv('list.csv', index=True)
#datetimes = np.array(datetimes)
for k,v in data.items():
data[k] = np.array(v)
return data
def plot(lats, lons, z = [], title = "Argo profiles south of 60S",
cbtitle = "Number of points in bin", vmax=None):
# setup north polar stereographic basemap.
# The longitude lon_0 is at 6-o'clock, and the
# latitude circle boundinglat is tangent to the edge
# of the map at lon_0. Default value of lat_ts
# (latitude of true scale) is pole.
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
m.drawcoastlines()
m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
plt.title("{} {}".format(len(lats), title))
x, y = m(lons, lats)
if len(z) == 0:
hh, locx, locy = np.histogram2d(x, y, bins=100)
# Sort the points by density, so that the densest points are plotted last
z = np.array([hh[np.argmax(a<=locx[1:]),np.argmax(b<=locy[1:])] for a,b in zip(x,y)])
#print(describe(z))
idx = z.argsort()
x, y, z = x[idx], y[idx], z[idx]
#m.imshow(heatmap, interpolation='bicubic', cmap="jet")
m.scatter(x, y, c=z, s=5, alpha=1, cmap="jet", vmax=vmax)
cb = plt.colorbar()
cb.set_label(cbtitle, rotation=270)
plt.show()
def plot_time(dts, title, label):
fig, ax = plt.subplots(1, 1, figsize=(20,10))
for i, dt in enumerate(dts):
#1 enumerate() method adds a counter to an iterable and returns it in a form of enumerate object.
# This enumerate object can then be used directly in for loops or be converted into a list of tuples using list() method.
k = label[i]
df = pd.DataFrame({k: pd.to_datetime(dt)})
print("df", df)
#1 Pandas was imported as pd
#2 Dataframe class - class pandas.DataFrame(data=None, index=None, columns=None, dtype=None, copy=False)
# Two-dimensional size-mutable, potentially heterogeneous tabular data structure with labeled axes (rows and columns).
# Arithmetic operations align on both row and column labels.
# Can be thought of as a dict-like container for Series objects
df = df.groupby(by=df[k].dt.date).count()
print ("df2", df)
#1 groupby - groups data by a certain specified variable type
df.plot(ax=ax, style='.', markersize=3)
ax.xaxis.set_major_locator(mdates.YearLocator())
ax.set_title(title, fontsize=18)
ax.set_xlabel("Time",fontsize=14)
ax.set_ylabel("Number of profiles",fontsize=14)
#ax.set_xlabel("")
ax.legend(fontsize=14)
plt.show()
argo = load_netcdf()
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:40: DeprecationWarning: tostring() is deprecated. Use tobytes() instead.
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:82: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
plot(argo["lat"], argo["lon"], vmax=100)
sst = [x[0] for x in argo["temperature"]]
plot(argo["lat"], argo["lon"], sst, title = "Argo SST", cbtitle = "Temperature °C")
max_depth = [np.nanmax(x) for x in argo["pressure"]]
plot(argo["lat"], argo["lon"], max_depth, title = "Argo Max Depth", cbtitle = "Depth (decibar pressure)")
/usr/local/lib/python3.6/dist-packages/numpy/core/_asarray.py:136: UserWarning: Warning: converting a masked element to nan. return array(a, dtype, copy=False, order=order, subok=True)
seal = load_netcdf("seal")
<ipython-input-389-b5e27fd2c44f>:19: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0 Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook` for f in tqdm(files):
df sorted
0 2004-02-05 22:38:00.000002
1 2004-02-06 02:49:59.999996
2 2004-02-06 06:08:59.999998
3 2004-02-06 13:57:00.000004
4 2004-02-06 23:25:00.000004
...
118762 2018-01-12 16:10:00.000001
118763 2018-01-12 22:30:00.000000
118764 2018-01-13 03:39:59.999998
118765 2018-01-13 22:10:00.000001
118766 2018-01-14 23:00:00.000003
Name: 0, Length: 118767, dtype: datetime64[ns]
plot(seal["lat"], seal["lon"], title = "Seal data south of 60S", vmax=100)
sst = [x[0] for x in seal["temperature"]]
plot(seal["lat"], seal["lon"], sst, title = "Seal SST", cbtitle = "Temperature °C")
max_depth = [np.nanmax(x) for x in seal["pressure"]]
plot(seal["lat"], seal["lon"], max_depth, title = "Seal Max Depth", cbtitle = "Depth (decibar pressure)")
print(min(seal["datetime"]), max(seal["datetime"]))
plot_time([seal["datetime"]], "Seal data dailystamps", ["Elephant seal"])
2004-02-05 22:38:00.000002 2018-01-14 23:00:00.000003 df Elephant seal 0 2007-03-03 03:24:00.000005 1 2007-03-03 03:42:59.999997 2 2007-03-03 06:00:59.999996 3 2007-03-03 17:40:59.999998 4 2007-03-03 17:51:59.999998 ... ... 118762 2007-09-07 09:19:59.999999 118763 2007-09-08 03:19:59.999999 118764 2007-09-08 03:19:59.999999 118765 2007-09-13 17:09:59.999998 118766 2007-09-13 17:09:59.999998 [118767 rows x 1 columns] df2 Elephant seal Elephant seal 2004-02-05 1 2004-02-06 4 2004-02-07 2 2004-02-08 2 2004-02-09 2 ... ... 2018-01-10 1 2018-01-11 1 2018-01-12 3 2018-01-13 2 2018-01-14 1 [4074 rows x 1 columns]
all_lats = np.concatenate((seal["lat"], argo["lat"]))
all_lons = np.concatenate((seal["lon"], argo["lon"]))
#plot(all_lats, all_lons, title="Argo + seal data south of 60S", vmax=100)
plot_time([seal["datetime"], argo["datetime"]], "Argo + seal data dailystamps", ["Seal", "Argo"])
df Seal 0 2007-03-03 03:24:00.000005 1 2007-03-03 03:42:59.999997 2 2007-03-03 06:00:59.999996 3 2007-03-03 17:40:59.999998 4 2007-03-03 17:51:59.999998 ... ... 118762 2007-09-07 09:19:59.999999 118763 2007-09-08 03:19:59.999999 118764 2007-09-08 03:19:59.999999 118765 2007-09-13 17:09:59.999998 118766 2007-09-13 17:09:59.999998 [118767 rows x 1 columns] df2 Seal Seal 2004-02-05 1 2004-02-06 4 2004-02-07 2 2004-02-08 2 2004-02-09 2 ... ... 2018-01-10 1 2018-01-11 1 2018-01-12 3 2018-01-13 2 2018-01-14 1 [4074 rows x 1 columns] df Argo 0 2010-08-24 22:56:17 1 2010-12-12 01:12:37 2 2010-12-21 23:00:22 3 2010-12-31 20:59:07 4 2011-01-10 19:09:02 ... ... 64373 2019-03-05 07:44:44 64374 2019-03-15 00:31:36 64375 2019-04-03 09:36:22 64376 2019-04-13 02:03:14 64377 2019-04-22 18:23:48 [64378 rows x 1 columns] df2 Argo Argo 2001-12-21 1 2002-01-04 1 2002-01-11 2 2002-01-18 2 2002-01-25 2 ... ... 2019-06-30 17 2019-07-01 11 2019-07-02 16 2019-07-03 5 2019-07-04 6 [5608 rows x 1 columns]
np.save("argo", argo)
np.save("seal", seal)
print(len(seal["lon"][:]), len(seal["pressure"][:]))
118767 118767
a3_sum = sum(argo["lat"]> -89.)
mask_out = argo["lat"] < -75
mask_in = (argo["lat"] < -60) & (argo["lat"] > -75)
a1_sum = sum(mask_out)
a2_sum = sum(mask_in)
print("Argo","\nA) Number of profiles to the south of 75S =",a1_sum,
"\nB) Number of profiles between 60 and 75S =", a2_sum,
"\nC) A+B =", (a2_sum+a1_sum),
"\nD) Total number of profiles to the south of 60S =", a3_sum)
s1_sum = sum(seal["lat"]> -89.)
mask_half_out = (seal["lat"] < -60) & (seal["lat"] > -89)
s2_sum = sum(mask_half_out)
mask_out = seal["lat"] < -75
mask_in = (seal["lat"] < -60) & (seal["lat"] > -75)
s3_sum = sum(mask_out)
s4_sum = sum(mask_in)
print("Seal","\nE) Number of profiles to the south of 75S =",s3_sum,
"\nF) Number of profiles between 60 and 75S =", s4_sum,
"\nG) E+F =", (s4_sum+s3_sum),
"\nH) Total number of profiles to the south of 60S =", s1_sum,
"\nI) Total number of profiles to the south of 60S =",s2_sum)
Argo A) Number of profiles to the south of 75S = 2571 B) Number of profiles between 60 and 75S = 61807 C) A+B = 64378 D) Total number of profiles to the south of 60S = 64378 Seal E) Number of profiles to the south of 75S = 3546 F) Number of profiles between 60 and 75S = 115221 G) E+F = 118767 H) Total number of profiles to the south of 60S = 118767 I) Total number of profiles to the south of 60S = 118767
# VERTICAL PROFILE SECTION (S,T,p)
from scipy.interpolate import interp1d
from scipy import interpolate
from scipy import signal
import pylab as plot
params = {'legend.fontsize': 14,
'legend.handlelength': 1,
'legend.labelspacing': 0.25,
'legend.handletextpad': 0.2,
'legend.frameon': False,
'legend.markerscale': 2.0,
'font.size': 20}
plot.rcParams.update(params)
def format_axes(fig):
for i, ax in enumerate(fig.axes):
ax.text(-1.5, 1950, "4.%d" % (i+1), va="center", ha="center", fontsize=12, color='black')
ax1.text(-1., 1950, "Total of Seal Profiles = {}".format(n_seal_profiles),
va="center", ha="left", fontsize=12, color='black')
ax3.text(-1., 1950, "Total of Argo Profiles = {}".format(n_argo_profiles),
va="center", ha="left", fontsize=12, color='black')
#ax1.text(2., 1750, 'Boundaries:\nLat (-{})-(-{})\nLon ({})-({})'.format(lat_i,lat_f,lon_i,lon_f),
# va="center", ha="center", fontsize=12, color='black')
ax.set_ylim(p_max, p_min);
ax.set_xlim(t_min, t_max)
ax.tick_params(axis="both", direction='inout', which='major', labelsize=16)
ax1.tick_params(labelbottom=False, labelleft=True)
ax2.tick_params(labelbottom=False, labelleft=False)
ax3.tick_params(labelbottom=True, labelleft=True)
ax4.tick_params(labelbottom=True, labelleft=False)
def format_axes_a(fig):
for i, ax in enumerate(fig.axes):
ax.text(-1.5, 1950, "6.%d" % (i+1), va="center", ha="center", fontsize=14, color='black')
ax.tick_params(axis="both", direction='inout', which='major', labelsize=16, labelbottom=True, labelleft=True)
ax.set_ylim(p_max, p_min);
ax.set_xlim(t_min, t_max)
ax1.set_ylabel("Depth (dbar)")
ax2.set_ylabel("Depth (dbar)") # , fontsize="20"
ax2.set_xlabel("Temperature (ºC)")
ax4.set_xlabel("Temperature (ºC)")
ax6.set_xlabel("Temperature (ºC)")
ax7.set_xlabel("Temperature (ºC)")
ax1.tick_params(labelbottom=False)
ax3.tick_params(labelbottom=False, labelleft=False)
ax4.tick_params(labelleft=False)
ax6.tick_params(labelleft=False)
ax5.tick_params(labelbottom=False,labelleft=False)
ax7.tick_params(labelbottom=True, labelleft=True)
def format_axes_b(fig):
for i, ax in enumerate(fig.axes):
ax.text(-1.5, 1950, "6.%d" % (i+1), va="center", ha="center", fontsize=14, color='black')
ax.tick_params(axis="both", direction='inout', which='major', labelsize=16, labelbottom=True, labelleft=True)
ax.set_ylim(p_max, p_min);
ax.set_xlim(t_min, t_max)
ax1.set_ylabel("Depth (dbar)", fontsize="20")
ax.set_xlabel("Temperature (ºC)")
#ax1.tick_params(labelbottom=True)
ax2.tick_params(labelleft=False)
ax3.tick_params(labelleft=False)
def format_axes_c(fig):
for i, ax in enumerate(fig.axes):
ax.text(-1.5, 1950, "6.%d" % (i+1), va="center", ha="center", fontsize=14, color='black')
ax.tick_params(axis="both", direction='inout', which='major', labelsize=16, labelbottom=True, labelleft=True)
ax.set_ylim(p_max, p_min);
ax.set_xlim(t_min, t_max)
ax1.set_ylabel("Depth (dbar)")
ax1.set_xlabel("Temperature (ºC)")
ax2.set_xlabel("Temperature (ºC)")
ax2.tick_params(labelbottom=True, labelleft=False)
def vert_prof_filter (vprofile, window, polyorder):
data = vprofile.interpolate(method='linear', limit= 200, limit_area="Inside")
#print("vertical profile ==> ", data)
#print("Vertical profile shape ==> ", data.shape)
y_smooth=signal.savgol_filter(data,
window_length=window, # window size used for filtering
polyorder=polyorder, # order of fitted polynomial
mode='nearest', cval=5), # Value to fill past the edges of the input if mode is ‘constant’. Default is 0.0.
#data_in_array = np.array(y_smooth)
##print("data_in_array.shape",data_in_array.shape)
##print("data_in_array ==> ", data_in_array)
return np.array(y_smooth)
# VERTICAL PROFILE SECTION (S,T,p)
# Defining Geographycal Boundaries
# -- 1º res
lat_n, lat_s = [65, 66] ### 0.5º res --> lat_n = 63.25 ; lat_s = 63.75 ; NO NEGATIVE
# n = northermost # s = southermost
lon_w, lon_e = [0.5, 1.5] # <-- 1º res ### 0.5º res --> #lon_w = 105.25 ; lon_e = 105.75;
# w = westermost # e = eastermost
# Grid system cell related to the geographycal boundaries
lat_grid, lon_grid = [6, 181]
#lat_n, lat_s = [63.5, 64.5]; lon_w, lon_e = [29.5, 30.5]; lat_grid, lon_grid = [4, 210]
fontsize=16
fontsize2=20
# defining the "area of interest"
mask_seal = ((seal["lat"] < -lat_n) & (seal["lat"] > -lat_s) & (seal["lon"] > lon_w) & (seal["lon"] < lon_e))
mask_argo = ((argo["lat"] < -lat_n) & (argo["lat"] > -lat_s) & (argo["lon"] > lon_w) & (argo["lon"] < lon_e))
# Maximum and minimum (S,T,p) in the plot
t_max, t_min = [3,-2]
s_max, s_min = [36., 33.]
p_max, p_min = [2000, 0]
depth = np.arange(p_min, p_max, 10)
xaxis = np.arange(t_min, t_max, 1)
print(xaxis)
# calculate the number of profiles within the "area of interest"
# In this case "np.flatnonzero" will return indices
# that are non-zero (or TRUE) in the flattened version of "mask_*"
n_seal_profiles = len(np.flatnonzero(mask_seal)) # Seal
n_argo_profiles = len(np.flatnonzero(mask_argo)) # Argo
#print("seal", len(np.flatnonzero(mask_seal)),np.flatnonzero(mask_seal))
#print("argo", len(np.flatnonzero(mask_argo)), np.flatnonzero(mask_argo))
print("Number of Seal Profiles in the window", len(np.flatnonzero(mask_seal)))
print("Number of Argo Profiles in the window", len(np.flatnonzero(mask_argo)))
print("seal", np.flatnonzero(mask_seal))
print("argo", np.flatnonzero(mask_argo))
[-2 -1 0 1 2] Number of Seal Profiles in the window 5 Number of Argo Profiles in the window 325 seal [25021 25023 25024 25025 25026] argo [ 9378 9379 9380 9381 9382 9383 9450 9451 9452 9453 9454 9455 9456 9457 9458 9459 9460 9461 9462 9463 9464 9465 9466 9467 9468 9469 9470 9471 9472 9473 9474 9475 9476 9477 9478 9479 9480 9481 9482 9483 9484 9485 9486 9487 9490 9491 9492 9493 9494 9495 9496 9497 9498 9499 9500 9501 9502 9503 9504 9505 9506 9507 9508 9509 9510 9511 9512 9513 9514 9515 9516 9517 9518 9519 9520 9521 9522 9523 9524 9525 9526 9527 9528 9529 9530 9531 9532 9533 9534 9535 9536 9537 9538 9539 9540 9541 9542 9543 9544 9545 9546 9547 9548 9549 9550 9551 9552 9553 9554 9555 9556 9557 9558 9559 9560 9561 9562 9563 9564 9565 9566 9567 9568 9569 9570 9571 9572 9573 9586 9587 9588 9589 9590 9591 9592 9593 9594 9595 9596 9597 9598 9599 9600 9601 9602 9603 9604 9605 9606 9607 9608 9609 9610 9611 9612 9613 9614 9615 9644 9645 9646 9647 9648 9649 9650 9651 9652 9653 9654 9655 9656 9657 9658 9659 9660 9661 9662 9663 9664 9665 9666 9667 9668 9669 9670 9671 9672 9673 9674 9675 9676 9677 9678 9679 9680 9681 9682 9683 9684 9685 9686 9687 9688 9689 9690 9691 17646 27215 27216 27223 27224 27225 27226 27227 27228 27229 27230 27231 27232 27233 27234 27235 27236 27237 27238 27239 27240 27241 27242 27243 27244 27245 27980 27981 27982 27983 27984 29317 29319 29320 29321 30596 30605 30606 32264 32265 32266 32267 32330 32331 32332 32333 32334 32335 32336 32337 32338 32339 32340 32341 32342 32343 32344 32345 32346 32347 32348 32349 32350 32351 32352 32353 32354 32355 32356 32357 32358 32359 32360 32361 32362 32363 32364 32365 32366 32367 32376 32377 32378 32379 32380 32381 32382 32383 32384 32385 32386 32387 32388 32389 32390 32391 32392 32393 37683 37684 37685 37686 37687 37688 37689 45646 45647 45648 45649 45650 45651 45652 51192 51193 51194 51204 51205 51206 51207] cabeca [False False False ... False False False]
for s in np.flatnonzero(mask_argo)[:10]:
print(s, argo["datetime"][s], argo["lat"][s], argo["lon"][s], argo["pressure"][s].mean(), argo["temperature"][s].mean())
#print(s,argomean[s],stdev[s])
9378 2012-01-02 10:58:52 -65.447 0.52 533.9277573529412 0.48338000727634806 9379 2012-01-02 10:58:52 -65.447 0.52 593.1396846064815 0.4775893599898727 9380 2012-01-09 11:01:21 -65.421 0.607 533.9314338235295 0.5161616306678921 9381 2012-01-09 11:01:21 -65.421 0.607 593.0285734953703 0.5081995151661061 9382 2012-01-16 11:16:39 -65.292 0.55 532.7249266144814 0.5112698848000244 9383 2012-01-16 11:16:39 -65.292 0.55 592.8871166087963 0.5051439779776113 9450 2012-09-10 19:25:09 -65.571 0.512 541.2585735586481 0.4319616359460425 9451 2012-09-10 19:25:09 -65.571 0.512 601.2847344483569 0.4317084657194469 9452 2012-09-17 18:51:43 -65.556 0.537 541.2632331013916 0.460987250326168 9453 2012-09-17 18:51:43 -65.556 0.537 601.2460754107981 0.4611112567740427
# VERTICAL PROFILE SECTION (S,T,p)
# Interpolate and plot "Individual" and "Mean" vertical profiles
#
# *** READ ME FIRST
# Firstly, must load "argo_temp_grid_withdepth" or "argo_sal_grid_withdepth"
fig = plt.figure(constrained_layout=True, figsize=(20,15))
gs = GridSpec(2, 4, figure=fig)
ax1 = fig.add_subplot(gs[0, 0]); ax3 = fig.add_subplot(gs[0, 1]); ax5 = fig.add_subplot(gs[0, 2]); ax7 = fig.add_subplot(gs[0, 3]);
ax2 = fig.add_subplot(gs[1, 0]); ax4 = fig.add_subplot(gs[1, 1]); ax6 = fig.add_subplot(gs[1, 2]); ax8 = fig.add_subplot(gs[1, 3]);
#ax7 = fig.add_subplot(gs[:, 3]); #ax7 = fig.add_subplot(gs[:, 3]) #fig.add_subplot(gs[:, 2:])
#ax7 = fig.add_subplot(gs[:, 3]) #fig.add_subplot(gs[:, 2:])
for i in np.flatnonzero(mask_seal):
t_ = seal["temperature"][i] #s_ = seal["salinity"][i]
p_ = seal["pressure"][i]
ax1.plot(t_, p_,'.', markersize=4)
ax1.legend(title='Seal profiles={}\nLat-{}-{}\nLon-{}-{}'.format(n_seal_profiles,lat_n, lat_s, lon_w, lon_e),
loc='lower right', title_fontsize=fontsize)
for i in np.flatnonzero(mask_argo):
t_ = argo["temperature"][i]
p_ = argo["pressure"][i]
ax2.plot(t_, p_,'.', markersize=4)
ax2.legend(title='Argo profiles={}\nLat-{}-{}\nLon-{}-{}'.format(n_argo_profiles,lat_n, lat_s, lon_w, lon_e),
loc='lower right', title_fontsize=fontsize)
#ax3.plot(seal_temp_grid_withdepth[lat_grid, lon_grid, :], depth, '.', markersize=6)
ax3.errorbar(seal_temp_grid_withdepth[lat_grid, lon_grid, :], depth,
xerr=seal_temp_grid_withdepth_stdev[lat_grid, lon_grid, :], fmt='.k', ecolor = 'blue', markersize=6, capsize=5)
ax3.legend(title='Seal mean\nLat ({})\nLon ({})'.format(lat_grid,lon_grid), loc='lower right', title_fontsize=fontsize)
#ax4.plot(argo_temp_grid_withdepth[lat_grid, lon_grid, :], depth, '.', markersize=6)
ax4.errorbar(argo_temp_grid_withdepth[lat_grid, lon_grid, :], depth,
xerr=argo_temp_grid_withdepth_stdev[lat_grid, lon_grid, :], fmt='.k', ecolor='blue', markersize=6, capsize=5)
ax4.legend(title='Argo mean\nLat ({})\nLon ({})'.format(lat_grid,lon_grid),
loc='lower right', title_fontsize=fontsize)
#interp_seal = interp_nans(seal_temp_grid_withdepth[:, lon_grid, :])
vprofile = pd.Series(seal_temp_grid_withdepth[lat_grid, lon_grid, :]) #print(vprofile)
window=9; polyorder=3;
fvprofile = vert_prof_filter(vprofile, window, polyorder)
ax5.plot(fvprofile[0,:], depth, '.', markersize=6)#, label="smooth w11 p.order 5")#, color='blue') #Savitzky-Golay
ax5.legend(title='Filtered mean\nLat ({})\nLon ({})'.format(lat_grid,lon_grid),
loc='lower right', title_fontsize=fontsize)
vprofile = pd.Series(argo_temp_grid_withdepth[lat_grid, lon_grid, :])
window=21; polyorder=5;
fvprofile = vert_prof_filter(vprofile, window, polyorder)
ax6.plot(fvprofile[0,:], depth, '.', markersize=6)#, label='Argo interpolated')
ax6.legend(title='Filtered mean\nLat ({})\nLon ({})'.format(lat_grid,lon_grid),
loc='lower right', title_fontsize=fontsize)
vprofile = pd.Series(combined_temp_grid_withdepth[lat_grid, lon_grid, :])
window=21; polyorder=5;
fvprofile = vert_prof_filter(vprofile, window, polyorder)
n_smooth, = ax7.plot(vprofile, depth, '.', markersize=4, color='blue',label='Combined mean')
v_smooth, = ax7.plot(fvprofile[0,:], depth, '.', color='green',markersize=6,
label='Filtered')
ax7.legend(handles=[n_smooth, v_smooth],loc='lower right', fontsize=fontsize)
ax8.plot(fvprofile[0,:], depth, '-k', markersize=10)
ax8.errorbar(combined_temp_grid_withdepth[lat_grid, lon_grid, :], depth,
xerr=combined_temp_grid_withdepth_stdev[lat_grid, lon_grid, :], fmt='.y', ecolor='blue', markersize=6, capsize=5)
ax8.legend(title='Combined mean\nLat ({})\nLon ({})'.format(lat_grid,lon_grid),
loc='lower right', title_fontsize=fontsize)
fig.suptitle("[Individual profiles] vs [mean profile] vs [Interpolated+smoothed profile]", fontsize=20)
format_axes_a(fig)
plt.show()
#print("dataframe",dataf)
No handles with labels found to put in legend. No handles with labels found to put in legend. No handles with labels found to put in legend. No handles with labels found to put in legend. No handles with labels found to put in legend. No handles with labels found to put in legend. No handles with labels found to put in legend.
# Lab Area for Vertical Profile Analysis (S,T,p)
from scipy.interpolate import interp1d
from scipy import interpolate
from scipy import signal
vprofile = pd.Series(combined_temp_grid_withdepth[lat_grid, lon_grid, :])
print("vprofile.shape",vprofile.shape)
y_smooth=signal.savgol_filter(vprofile,
window_length=11, # window size used for filtering
polyorder=5, # order of fitted polynomial
mode='interp', # ‘mirror’, ‘constant’, ‘nearest’, ‘wrap’ or ‘interp’. This determines the type of extension to use for the padded signal to which the filter is applied.
cval=5), # Value to fill past the edges of the input if mode is ‘constant’. Default is 0.0.
#print("vprofile",vprofile,"y_smooth", y_smooth)
#print("y_smooth", y_smooth);print("length y_smooth", len(y_smooth));
data_in_array = np.array(y_smooth)
#print("data_in_array", data_in_array);print("data_in_array.shape", data_in_array.shape)
ss = pd.Series(y_smooth) # it is creating a serie (or list) of lists
#print ("ss --> ",ss)
fig = plt.figure(figsize=(10,20))
ax = fig.add_subplot(111)
ax.plot(combined_temp_grid_withdepth[lat_grid, lon_grid, :], depth,
'--', color='green', markersize=4, label="Combined Mean Field")
ax.plot(data_in_array[0,:], depth, '--',
label="Smooth Window-11 poly.order.5", color='blue') #Savitzky-Golay
ax.errorbar(combined_temp_grid_withdepth[lat_grid, lon_grid, :], depth,
xerr=combined_temp_grid_withdepth_stdev[lat_grid, lon_grid, :], fmt='.y', ecolor='blue', markersize=6, capsize=5)
format_axes_b(fig)
plt.show()
vprofile.shape (200,)
# Lab Area for Vertical Profile Analysis (S,T,p)
# Confront data vs filtered data
from scipy.interpolate import interp1d
from scipy import interpolate
from scipy import signal
f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, sharex=True, figsize=(12,8))
for i in np.flatnonzero(mask_argo):
temp_argo = argo["temperature"][i]
pressure_argo = argo["pressure"][i]
y_smooth=signal.savgol_filter(temp_argo,
window_length=21, # window size used for filtering
polyorder=5, # order of fitted polynomial
mode='nearest',
cval=5), # Value to fill past the edges of the input if mode is ‘constant’. Default is 0.0.
data_in_array = np.array(y_smooth)
v_smooth, = ax2.plot(data_in_array[0,:], pressure_argo, '.', markersize=2,
label='Filtered\nProfiles') #Savitzky-Golay
n_smooth, = ax1.plot(temp_argo, pressure_argo,'.', markersize=2,
label='Total of\nProfiles {}'.format(len(np.flatnonzero(mask_argo))))
ax1.errorbar(argo_temp_grid_withdepth[lat_grid, lon_grid, :], depth,
xerr=argo_temp_grid_withdepth_stdev[lat_grid, lon_grid, :],
fmt='.k', ecolor='black', markersize=6, capsize=3)
ax2.legend(handles=[v_smooth], loc='upper right',fontsize=12)
ax1.legend(handles=[n_smooth],loc='upper right',fontsize=12)
#ax1.legend(title='Argo profiles={}\nLat-{}-{}\nLon-{}-{}'.format(n_argo_profiles,lat_n, lat_s, lon_w, lon_e),
# loc='upper right', title_fontsize=fontsize)
format_axes_c(f)
plt.show()
for f in glob.glob("data/ctd/**/*.nc", recursive=True):
d = Dataset(f)
try:
temp = d.variables["Temperature"][:]
if np.nanmin(temp) < -1000:
print(f, np.nanmin(temp), np.nanmean(temp), np.nanmax(temp), temp[temp < -1000])
except:
pass
data/ctd/wod_019543360O.nc -10000000000.0 -33660590.0 1.29 [-1.e+10 -1.e+10 -1.e+10 -1.e+10 -1.e+10 -1.e+10 -1.e+10 -1.e+10 -1.e+10 -1.e+10 -1.e+10 -1.e+10] data/ctd/wod_017399145O.nc -10000000000.0 -3841720.0 2.0493 [-1.e+10] data/ctd/wod_013620747O.nc -10000000000.0 -10095910.0 1.9588 [-1.e+10 -1.e+10] data/ctd/wod_017730902O.nc -10000000000.0 -75187970.0 1.0537 [-1.e+10] data/ctd/wod_013620743O.nc -10000000000.0 -3040437.5 1.4328 [-1.e+10] data/ctd/wod_019543407O.nc -10000000000.0 -5149330.5 1.704 [-1.e+10 -1.e+10] data/ctd/wod_016492374O.nc -10000000000.0 -2214839.8 0.5483 [-1.e+10] data/ctd/wod_016492449O.nc -10000000000.0 -25380710.0 2.7721 [-1.e+10] data/ctd/wod_015561197O.nc -10000000000.0 -3935457.2 2.0305 [-1.e+10]
def load_netcdf_grid(path = "Argo_South_60", resolution = 1, target_var = "temp",
with_depth = False, filter_month = None, filter_season = None,
filter_start = datetime(2000,1,1), filter_end = datetime(2021,1,1), mean = None):
if with_depth:
grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
gridcounts = np.zeros(shape=[15 * resolution, 360 * resolution, 200])
else:
grid = np.full(shape=[15 * resolution, 360 * resolution], fill_value=np.nan)
files = glob.glob("data/{}/**/*.nc".format(path), recursive=True)
#1 glob.glob - returns a list of all the files in a given folder with the extension 'nc'
#2 recursive=true - keeps on looking for all the files till the end, or the last file that exist
# creating lists
filename2 = []
nbprof2 = []
filename2_maskout = []
## End of "creating lists"
for f in tqdm(files[:]):
# tqdm - make your loops show a progress meter
d = Dataset(f)
#1 DATASET - Creating/Opening/Closing a netCDF file.
# This is the method used to open an existing netCDF file.
lat = d.variables["LATITUDE"][:]
#mask = (lat < -60) & (lat > -74.75)
mask = (lat < -60) & (lat > -74.5)
# Counting number of profiles per NetCDF files
nbprof = [len(lat)] # number of profiles per file
nbprof2.extend(nbprof)
# filename = filename
filename = ["Filename {}".format(f)]
filename2.extend(filename)
#Create a DataFrame: filename + number of profiles per file
profiles_per_files = {'Filename': filename2, 'Number of profiles per file': nbprof2}
new_df = pd.DataFrame(data=profiles_per_files)
#Save "new_df" as csv file:
new_df.to_csv('filename_nbprof.csv', index=False)
## End of "Counting number of profiles per NetCDF files"
juld = d.variables["JULD"][:]
units = d.variables["JULD"].getncattr('units')
dates = num2date(juld, units, "standard")
datemask = (dates > filter_start) & (dates < filter_end)
mask &= datemask
if filter_month:
datemask = np.array([d.month == filter_month for d in dates])
mask &= datemask
if filter_season:
if filter_season == "Summer":
datemask = np.array([d.month >= 10 or d.month <= 3 for d in dates])
elif filter_season == "Winter":
datemask = np.array([d.month > 3 and d.month < 10 for d in dates])
mask &= datemask
if any(mask):
lat = np.round((np.abs(lat[mask]) - 60) * resolution).astype(int)
lon = d.variables["LONGITUDE"][mask]
lon = np.round((lon + 180) * resolution).astype(int)
lon[lon == 360 * resolution] = 0
if with_depth:
#pres = np.round(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
pres = np.floor(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
temp = d.variables["TEMP_ADJUSTED"][mask]
if target_var == "sal":
try:
sal = d.variables["PSAL_ADJUSTED"][mask]
except:
print("No salinity for {}".format(f))
continue
else:
pres = d.variables["PRES_ADJUSTED"][mask]
temp = d.variables["TEMP_ADJUSTED"][mask, 0]
for x in np.unique(lon):
for y in np.unique(lat):
if with_depth:
ptmask = (lon == x) & (lat == y)
for z in np.unique(pres[ptmask]):
if z>= 200 or np.ma.is_masked(z):
continue
depthmask = pres[ptmask] == z
if target_var == "temp":
values_at_pt = temp[ptmask][depthmask].compressed()
elif target_var == "sal":
values_at_pt = sal[ptmask][depthmask].compressed()
if mean is not None:
mean_at_pt = mean[y, x, z]
grid[y, x, z] = np.nansum((grid[y, x, z],np.nansum(np.square(values_at_pt - mean_at_pt))))
else:
grid[y, x, z] = np.nansum((grid[y, x, z], np.nansum(values_at_pt)))
gridcounts[y, x, z] += len(values_at_pt)
else:
ptmask = (lon == x) & (lat == y)
if target_var == "n_point":
v = np.sum(ptmask)
if not np.ma.is_masked(v):
if v != 0:
#1 added a if statement checking if number of points is greater than 0
grid[y, x] = np.nansum((grid[y, x], v))
elif target_var == "temp":
temps_at_pt = temp[ptmask]
v = np.ma.mean(temps_at_pt)
# 1 MEAN - numpy.ma.mean(self, axis=None, dtype=None, out=None, keepdims=<no value>)
# Returns the average of the array elements along given axis.
# Masked entries are ignored, and result elements which are not finite
# will be masked.
sd = np.ma.std(temps_at_pt)
# numpy.std -->> Standard deviation
if not np.ma.is_masked(v):
grid[y, x] = np.nanmean((grid[y, x], v))
elif target_var == "depth":
pres_at_pt = pres[ptmask]
if len(pres_at_pt):
v = np.ma.max(pres[ptmask])
if not np.ma.is_masked(v):
grid[y, x] = np.nanmax((grid[y, x], v))
if with_depth:
grid /= gridcounts
print(np.max(gridcounts))
print(np.nanmin(grid), np.nanmean(grid), np.nanmax(grid))
if mean is not None:
grid = np.sqrt(grid)
return grid
def load_ctd_grid(resolution = 1, target_var = "temp",
with_depth = False, filter_month = None, filter_season = None,
filter_start = datetime(2000,1,1), filter_end = datetime(2021,1,1), mean = None):
if with_depth:
grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
gridcounts = np.zeros(shape=[15 * resolution, 360 * resolution, 200])
else:
grid = np.full(shape=[15 * resolution, 360 * resolution], fill_value=np.nan)
files = glob.glob("data/ctd/**/*.nc", recursive=True)
#1 glob.glob - returns a list of all the files in a given folder with the extension 'nc'
#2 recursive=true - keeps on looking for all the files till the end, or the last file that exist
for f in tqdm(files[:]):
# tqdm - make your loops show a progress meter
d = Dataset(f)
#1 DATASET - Creating/Opening/Closing a netCDF file.
# This is the method used to open an existing netCDF file.
lat = d.variables["lat"][:]
if lat.ndim > 0:
continue
if lat > -60 or lat < -74.5:
continue
time = d.variables["time"]
units = time.getncattr('units')
date = num2date(time[:], units, "standard")
if date < filter_start or date > filter_end:
continue
if filter_month:
if date.month != filter_month:
continue
if filter_season:
if filter_season == "Summer":
if not (date.month >= 10 or date.month <= 3):
continue
elif filter_season == "Winter":
if not (date.month > 3 and date.month < 10):
continue
lat = np.round((np.abs(lat) - 60) * resolution).astype(int)
lon = d.variables["lon"][:]
lon = np.round((lon + 180) * resolution).astype(int)
if lon == (360 * resolution):
lon = 0
x = lon
y = lat
temp = d.variables["Temperature"][:]
temp[temp == -1e10] = np.nan # Way too cold!
if with_depth:
#pres = np.round(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
pres = np.floor(d.variables["Pressure"] / 10).astype(int)
if target_var == "sal":
try:
sal = d.variables["PSAL_ADJUSTED"]
except:
print("No salinity for {}".format(f))
continue
else:
pres = d.variables["Pressure"]
temp = temp[0]
if with_depth:
for z in np.unique(pres):
if z>= 200 or np.ma.is_masked(z):
continue
depthmask = pres == z
if target_var == "temp":
values_at_pt = temp[depthmask].compressed()
elif target_var == "sal":
values_at_pt = sal[depthmask].compressed()
if mean is not None:
mean_at_pt = mean[y, x, z]
grid[y, x, z] = np.nansum((grid[y, x, z],np.nansum(np.square(values_at_pt - mean_at_pt))))
else:
grid[y, x, z] = np.nansum((grid[y, x, z], np.nansum(values_at_pt)))
gridcounts[y, x, z] += len(values_at_pt)
else:
if target_var == "n_point":
grid[y, x] = np.nansum((grid[y, x], 1))
elif target_var == "temp":
temps_at_pt = temp
v = np.ma.mean(temps_at_pt)
# 1 MEAN - numpy.ma.mean(self, axis=None, dtype=None, out=None, keepdims=<no value>)
# Returns the average of the array elements along given axis.
# Masked entries are ignored, and result elements which are not finite
# will be masked.
sd = np.ma.std(temps_at_pt)
# numpy.std -->> Standard deviation
if not np.ma.is_masked(v):
grid[y, x] = np.nanmean((grid[y, x], v))
elif target_var == "depth":
pres_at_pt = pres
if len(pres_at_pt):
v = np.ma.max(pres[ptmask])
if not np.ma.is_masked(v):
grid[y, x] = np.nanmax((grid[y, x], v))
if with_depth:
grid /= gridcounts
print(np.max(gridcounts))
print(np.nanmin(grid), np.nanmean(grid), np.nanmax(grid))
if mean is not None:
grid = np.sqrt(grid)
return grid
# This version seems a lot more inefficient for some reason
def load_netcdf_grid_new(platform = "argo", resolution = 1, target_var = "temp", filter_month = None,
filter_season = None, filter_start = datetime(2000,1,1), filter_end = datetime(2020,1,1)):
grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
if platform == "argo":
data = argo
elif platform == "seal":
data = seal
elif platform == "both":
data = {}
for k in argo:
data[k] = argo[k] + seal[k]
mask = (data["lat"] < -60) & (data["lat"] > -74.5)
datemask = (data["datetime"] > filter_start) & (data["datetime"] < filter_end)
mask &= datemask
if filter_month:
datemask = np.array([d.month == filter_month for d in data["datetime"]])
mask &= datemask
if filter_season:
if filter_season == "Summer":
datemask = np.array([d.month >= 10 or d.month <= 3 for d in data["datetime"]])
elif filter_season == "Winter":
datemask = np.array([d.month > 3 and d.month < 10 for d in data["datetime"]])
mask &= datemask
if any(mask):
lat = np.round((np.abs(data["lat"][mask]) - 60) * resolution).astype(int)
lon = data["lon"][mask]
lon = np.round((lon + 180) * resolution).astype(int)
lon[lon == 360 * resolution] = 0
depth = np.floor(data["pressure"][mask] / 10).compressed().astype(int)
grid[lat, lon, depth] = np.ma.mean(data["temperature"][mask])
return grid
def plot_grid(grid, resolution = 1, title="Ocean data", cmap="nipy_spectral", cbtitle="Temperature °C", vmin=None, vmax=None):
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
m.drawcoastlines()
m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, -80, -5))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
#plt.title("{}".format(title, 1/resolution),fontsize = 16)
#plt.title("\n{} at {}° resolution".format(title, 1/resolution),fontsize = 16)
x = np.arange(-180, 180.1, 1/resolution)
y = np.arange(-60, -75.1, 1/-resolution)
x, y = np.meshgrid(x, y)
x, y = m(x, y)
if cbtitle=="Number of profiles":
pcm = m.pcolormesh(x, y, grid, cmap=cmap, norm=colors.LogNorm(vmin, vmax))
# #norm=colors.PowerNorm(gamma=0.25))
# print("vmin",vmin,"vmax",vmax)
#ax.xaxis.grid(True, which='major')
cb = plt.colorbar(pcm, extend='max',shrink=0.785, pad=0.05 , orientation="horizontal")
#cb.locator=ticker.LogLocator(base=10)
cb.set_label(cbtitle, labelpad= 10, rotation=0, fontsize = 14)
else:
m.pcolormesh(x, y, grid, cmap=cmap, vmin=vmin, vmax=vmax) #"colorcet" "nipy_spectral" "ocean" "gist_ncar" "cubehelix" "jet" "PuOr" "RdBu" "rainbow"
#m.pcolormesh(x, y, grid, cmap=cmap, vmin=vmin, vmax=vmax) #"colorcet" "nipy_spectral" "ocean" "gist_ncar" "cubehelix" "jet" "PuOr" "RdBu" "rainbow"
cb = plt.colorbar(shrink=0.785, pad=0.05 , orientation="horizontal")
cb.set_label(cbtitle, labelpad= 10, rotation=0, fontsize = 14)
#m.pcolormesh(x, y, grid, cmap="jet", vmin=vmin, vmax=vmax)
plt.show()
def interp_nans(grid, method='linear'):
x = np.arange(0, grid.shape[1])
y = np.arange(0, grid.shape[0])
mask = np.isnan(grid)
xx, yy = np.meshgrid(x, y)
x1 = xx[~mask]
y1 = yy[~mask]
newarr = grid[~mask]
interp = interpolate.griddata((x1, y1), newarr.ravel(), (xx, yy), method=method)
if method == "linear":
interp = interp_nans(interp, "nearest")
return interp
ctd_grid = load_ctd_grid()
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:153: DeprecationWarning: tostring() is deprecated. Use tobytes() instead. /usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:217: RuntimeWarning: Mean of empty slice
print(np.nanmin(ctd_grid), np.nanmean(ctd_grid), np.nanmax(ctd_grid))
plot_grid(ctd_grid)
-1.8797187823802233 -0.21100522311239855 5.473199844360352
# Titles and labels for figures
a_title, s_title = ["Argo Data", "Seal Data"]
comb_title = "Argo+Seal Data"
interp_title = "Interpolated Argo+Seal Data"
smf = "Salinity mean field "
tmf = "Temperature mean field "
dmf = "Density mean field "
param_t, param_s = ["Temperature", "Salinity"]
param_d = "Density"
# Maximum and Minimum - temperature and Salinity
smax, smin = [34.8, 32.5]
tmax, tmin = [7, -2.5]
rhomin, rhomax = [1028, 1026.5]
# Grid Resolution
resolution = 1
argo_pt_grid = load_netcdf_grid(target_var="n_point", with_depth = False, resolution = resolution)
seal_pt_grid = load_netcdf_grid(path="seal", target_var="n_point", with_depth = False, resolution = resolution)
<ipython-input-283-7a4d97c2e055>:20: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0 Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook` for f in tqdm(files[:]):
combined_pt_grid = np.nansum((argo_pt_grid, seal_pt_grid), axis=0)
combined_pt_grid[np.isnan(argo_pt_grid) & np.isnan(seal_pt_grid)] = np.nan
print(argo_pt_grid.shape,seal_pt_grid.shape,combined_pt_grid.shape)
(15, 360) (15, 360) (15, 360)
# Lab log numbers
xx2 = np.nanmin(seal_pt_grid),np.nanmax(seal_pt_grid); print("xx2",xx2)
xx4 = np.log2(xx2); print("xx4",xx4)
xx6 = np.floor(xx4); print("xx6",xx6)
#xx7 = np.ceil(xx6); print(xx7)
lev_exp = np.arange(np.floor(np.log2(np.nanmin(seal_pt_grid))),
np.ceil(np.log2(np.nanmax(seal_pt_grid))+1))
print("lev_exp",lev_exp)
levs = np.power(2, lev_exp)
print("levs=", levs)
print("levs min and max",levs.min(),levs.max())
tickList = []
for i in range(0,len(lev_exp)):
value = float(levs[i])
tickList.append(value)
print("Float Levs [i]",tickList)
xx2 (1.0, 6015.0) xx4 [ 0. 12.55434902] xx6 [ 0. 12.] lev_exp [ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.] levs= [1.000e+00 2.000e+00 4.000e+00 8.000e+00 1.600e+01 3.200e+01 6.400e+01 1.280e+02 2.560e+02 5.120e+02 1.024e+03 2.048e+03 4.096e+03 8.192e+03] levs min and max 1.0 8192.0 Float Levs [i] [1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0, 512.0, 1024.0, 2048.0, 4096.0, 8192.0]
#plot_grid(argo_pt_grid, title=a_title + "\nNumber of vertical profiles", cmap="jet", cbtitle="Number of profiles", vmin=1,vmax=80, resolution = resolution)
#plot_grid(seal_pt_grid, title=s_title + "\nNumber of vertical profiles", cmap="jet", cbtitle="Number of profiles", vmin=1,vmax=100, resolution = resolution)
#plot_grid(combined_pt_grid, title=comb_title + "\nNumber of vertical profiles", cmap="jet", cbtitle="Number of profiles", vmin=1,vmax=100, resolution = resolution)
import matplotlib.colors as colors
import matplotlib.cbook as cbook
from matplotlib.colors import SymLogNorm, LogNorm
plot_grid(argo_pt_grid, cmap="jet",
cbtitle="Number of profiles", vmin=np.floor(np.nanmin(argo_pt_grid)),
vmax=np.ceil(np.nanmax(argo_pt_grid)),resolution = resolution)
plot_grid(seal_pt_grid, cmap="jet",
cbtitle="Number of profiles", vmin=np.floor(np.nanmin(seal_pt_grid)),
vmax=np.nanmax(seal_pt_grid), resolution = resolution)
c:\users\fvie285\appdata\local\programs\python\python38\lib\site-packages\matplotlib\colors.py:1171: RuntimeWarning: invalid value encountered in less_equal mask |= resdat <= 0
c:\users\fvie285\appdata\local\programs\python\python38\lib\site-packages\matplotlib\colors.py:1171: RuntimeWarning: invalid value encountered in less_equal mask |= resdat <= 0
argo_temp_grid = load_netcdf_grid(resolution = resolution)
seal_temp_grid = load_netcdf_grid("seal", resolution = resolution)
print(argo_temp_grid.shape,seal_temp_grid.shape)
<ipython-input-254-7a4d97c2e055>:20: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0 Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook` for f in tqdm(files[:]):
c:\users\fvie285\appdata\local\programs\python\python38\lib\site-packages\numpy\ma\core.py:5215: RuntimeWarning: Mean of empty slice. result = super(MaskedArray, self).mean(axis=axis, c:\users\fvie285\appdata\local\programs\python\python38\lib\site-packages\numpy\core\_methods.py:153: RuntimeWarning: invalid value encountered in true_divide ret = um.true_divide( c:\users\fvie285\appdata\local\programs\python\python38\lib\site-packages\numpy\ma\core.py:5293: RuntimeWarning: Degrees of freedom <= 0 for slice ret = super(MaskedArray, self).var(axis=axis, dtype=dtype, out=out, c:\users\fvie285\appdata\local\programs\python\python38\lib\site-packages\numpy\core\_methods.py:185: RuntimeWarning: invalid value encountered in true_divide arrmean = um.true_divide( c:\users\fvie285\appdata\local\programs\python\python38\lib\site-packages\numpy\core\_methods.py:206: RuntimeWarning: invalid value encountered in true_divide ret = um.true_divide( <ipython-input-254-7a4d97c2e055>:110: RuntimeWarning: Mean of empty slice grid[y, x] = np.nanmean((grid[y, x], v))
(15, 360) (15, 360)
combined_temp_grid = np.nanmean((argo_temp_grid, seal_temp_grid), axis=0)
print(combined_temp_grid.shape)
(15, 360)
<ipython-input-258-e15be8551138>:1: RuntimeWarning: Mean of empty slice combined_temp_grid = np.nanmean((argo_temp_grid, seal_temp_grid), axis=0)
#title = "Seal data gridded at {}\n0-10 decibar mean temperature\nMean across all points = {}".format(1/resolution,sealmean)
title = a_title + "\n0-10 decibars" + mtf
plot_grid(argo_temp_grid, title=title, resolution = resolution)
title = s_title + "\n0-10 decibar" + mtf
plot_grid(seal_temp_grid, title=title, resolution = resolution)
title = comb_title + "\n0-10 decibar" + mtf
plot_grid(combined_temp_grid, title=title, resolution = resolution)
interp = interp_nans(combined_temp_grid)
#title = "Argo+Seal data gridded at {}\n0-10 decibar mean temperature\nMean across all points = {}".format(1/resolution,interpmean)
title = interp_title + "\n0-10 decibar" + mtf
plot_grid(interp, resolution = resolution, title=title)
argo_max_depth = load_netcdf_grid(target_var = "depth", resolution = resolution)
seal_max_depth = load_netcdf_grid("seal", target_var = "depth", resolution = resolution)
#print(argo_max_depth.shape,seal_max_depth.shape)
plot_grid(argo_max_depth, resolution = resolution, title="Argo data Max depth\n", cbtitle="Depth (dbar)")
plot_grid(seal_max_depth, resolution = resolution, title="Seal data Max depth\n", cbtitle="Depth (dbar)", vmax=2000)
# Monthly temperature field
for month in range(1, 13):
try:
argo_temp = np.load(f"argo_temp_grid_{month}.npy")
seal_temp = np.load(f"seal_temp_grid_{month}.npy")
combined_temp = np.load(f"combined_temp_grid_{month}.npy")
except:
argo_temp = load_netcdf_grid(filter_month = month, with_depth=True)
seal_temp = load_netcdf_grid("seal", filter_month = month, with_depth=True)
combined_temp = np.nanmean((argo_temp, seal_temp), axis=0)
np.save("argo_temp_grid_{}".format(month), argo_temp)
np.save("seal_temp_grid_{}".format(month), seal_temp)
np.save("combined_temp_grid_{}".format(month), combined_temp)
interp_temp = interp_nans(combined_temp[:,:,0])
argomean = np.round(np.nanmean(argo_temp), 2)
sealmean = np.round(np.nanmean(seal_temp), 2)
combmean = np.round(np.nanmean(combined_temp), 2)
interpmean = np.round(np.nanmean(interp_temp), 2)
# #title = a_title + "\n0-10 decibar mean temperature for {}\nMean across all points = {}".format(calendar.month_abbr[month], argomean)
title = a_title + "\n0-10 decibars mean temperature for {}".format(calendar.month_abbr[month])
#plot_grid(argo_temp[:,:,0], title=title)
# #title = "Seal data\n0-10 decibar mean temperature for {}\nMean across all points = {}".format(calendar.month_abbr[month], sealmean)
title = s_title + "\n0-10 decibar mean temperature for {}".format(calendar.month_abbr[month])
#plot_grid(seal_temp[:,:,0], title=title)
title = comb_title + "\n0-10 decibar mean temperature for {}".format(calendar.month_abbr[month])
# title = "Argo+Seal data\n0-10 decibar mean temperature for {}\nMean across all points = {}".format(calendar.month_abbr[month], combmean)
plot_grid(combined_temp[:,:,0], title=title)
title = interp_title + "\n0-10 decibar mean for {}".format(calendar.month_abbr[month])
#plot_grid(interp_temp, title=title)
# Load / Create the Lat-Long-Depth Grid of the Seasonal Mean Temperature Field
# Plot Horizontal Sections of the Seasonal Mean Temperature Field
resolution = 1
#dbars = [0, 5, 10, 15, 20, 30, 50, 75, 100, 125, 150, 175, 190, 199]
title="Temperature (ºC)"
dbars = [0]
tmin, tmax = [-2.5, 6];
for season in ["Winter", "Summer"]:
try:
#argo_pts_grid = np.load(f"argo_pts_grid_{season}.npy")
argo_temp = np.load(f"argo_temp_grid_withdepth_{season}.npy")
seal_temp = np.load(f"seal_temp_grid_withdepth_{season}.npy")
combined_temp = np.load(f"combined_temp_grid_withdepth_{season}.npy")
except:
argo_temp = load_netcdf_grid(filter_season = season, with_depth=True, resolution = resolution)
seal_temp = load_netcdf_grid("seal", filter_season = season, with_depth=True, resolution = resolution)
combined_temp = np.nanmean((argo_temp, seal_temp), axis=0)
np.save("argo_temp_grid_withdepth_{}".format(season), argo_temp)
np.save("seal_temp_grid_withdepth_{}".format(season), seal_temp)
np.save("combined_temp_grid_withdepth_{}".format(season), combined_temp)
interp_temp = interp_nans(combined_temp[:,:,0])
argomean = np.round(np.nanmean(argo_temp), 2)
sealmean = np.round(np.nanmean(seal_temp), 2)
combmean = np.round(np.nanmean(combined_temp), 2)
interpmean = np.round(np.nanmean(interp_temp), 2)
for db in dbars:
print("valor minimo e maximo", np.nanmin(combined_temp[:,:,db]),np.nanmax(combined_temp[:,:,db]))
#dbs = "between {}-{} dbar ".format(db*10, (db+1)*10)
#title = a_title + "\n" + mtf + dbs + "for {}".format(season)
#plot_grid(argo_temp[:,:,db], title=title, resolution = resolution, vmin=tmin, vmax=tmax)
#title = "Argo data 0-10 decibar mean temperature for {}\nMean across all points = {}".format(season, argomean)
#title = s_title + "\n" + mtf + dbs + "for {}".format(season)
#plot_grid(seal_temp[:,:,db], title=title, resolution = resolution, vmin=tmin, vmax=tmax)
#title = comb_title + "\n" + mtf + dbs + "for {}".format(season)
plot_grid(combined_temp[:,:,db], title=title, resolution = resolution, vmin=tmin, vmax=tmax)
# Load / Create the Lat-Long-Depth Grid for the SEASONAL Mean SALINITY FIELD
# Plot Horizontal Sections of the Seasonal Mean Salinity Field
resolution = 1
#dbars = [0, 5, 10]
smin, smax = [32.5, 35]
dbars = [0]
title="Salinity"
for season in ["Winter","Summer"]:
try:
#argo_pts_grid = np.load(f"argo_pts_grid_{season}.npy")
argo_sal = np.load(f"argo_sal_grid_withdepth_{season}.npy")
seal_sal = np.load(f"seal_sal_grid_withdepth_{season}.npy")
combined_sal = np.load(f"combined_sal_grid_withdepth_{season}.npy")
except:
argo_sal = load_netcdf_grid(target_var = "sal", filter_season = season, with_depth=True, resolution = resolution)
seal_sal = load_netcdf_grid("seal", target_var = "sal", filter_season = season, with_depth=True, resolution = resolution)
combined_sal = np.nanmean((argo_sal, seal_sal), axis=0)
np.save("argo_sal_grid_withdepth_{}".format(season), argo_sal)
np.save("seal_sal_grid_withdepth_{}".format(season), seal_sal)
np.save("combined_sal_grid_withdepth_{}".format(season), combined_sal)
interp_sal = interp_nans(combined_sal[:,:,0])
argo_s_mean = np.round(np.nanmean(argo_sal), 2)
seal_s_mean = np.round(np.nanmean(seal_sal), 2)
comb_s_mean = np.round(np.nanmean(combined_sal), 2)
interp_s_mean = np.round(np.nanmean(interp_sal), 2)
for db in dbars:
dbs = "between {}-{} dbar ".format(db*10, (db+1)*10)
#title = a_title + "\n" + msf + dbs + "for {}".format(season)
print("valor minimo e maximo", np.nanmin(combined_sal[:,:,db]),np.nanmax(combined_sal[:,:,db]))
#plot_grid(argo_sal[:,:,db], title=title, cbtitle="Salinity", resolution = resolution, vmin=smin, vmax=smax)
#title = "Argo data 0-10 decibar mean temperature for {}\nMean across all points = {}".format(season, argomean)
#title = s_title + "\n" + msf + dbs + "for {}".format(season)
#plot_grid(seal_sal[:,:,db], title=title, cbtitle="Salinity", resolution = resolution, vmin=smin, vmax=smax)
#title = comb_title + "\n" + msf + dbs + "for {}".format(season)
plot_grid(combined_sal[:,:,db], title=title, cbtitle="Salinity", resolution = resolution, vmin=smin, vmax=smax)
print(argo_temp_grid_withdepth.shape, seal_temp_grid_withdepth.shape, argo_temp_grid_withdepth.shape)
#print(max(argo_temp[:,:,0]))
#print(argo_lats.shape, argo_lons.shape)
#print(seal_lats.shape, seal_lons.shape)
#print(argo_lats)
# Load / Create the Lat-Long-Depth Grid of the Climatological Temperature Field
resolution = 1
try:
argo_temp_grid_withdepth = np.load("argo_temp_grid_withdepth.npy")
seal_temp_grid_withdepth = np.load("seal_temp_grid_withdepth.npy")
combined_temp_grid_withdepth = np.load("combined_temp_grid_withdepth.npy")
except:
argo_temp_grid_withdepth = load_netcdf_grid(with_depth=True, resolution = resolution)
seal_temp_grid_withdepth = load_netcdf_grid("seal", with_depth=True, resolution = resolution)
combined_temp_grid_withdepth = load_netcdf_grid("*", with_depth=True, resolution = resolution)
np.save("argo_temp_grid_withdepth", argo_temp_grid_withdepth)
np.save("seal_temp_grid_withdepth", seal_temp_grid_withdepth)
np.save("combined_temp_grid_withdepth", combined_temp_grid_withdepth)
<ipython-input-317-7a4d97c2e055>:20: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0 Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook` for f in tqdm(files[:]):
1098.0 -1.9201643807547433 0.9862097505275941 5.1874998807907104
9537.0 -1.9967581033706665 0.38976182178951635 5.692234516143799
<ipython-input-317-7a4d97c2e055>:119: RuntimeWarning: invalid value encountered in true_divide grid /= gridcounts
9537.0 -1.9967581033706665 0.9278911659997323 5.1874998807907104
resolution = 1
try:
argo_temp_grid_withdepth_stdev = np.load("argo_temp_grid_withdepth_stdev.npy")
seal_temp_grid_withdepth_stdev = np.load("seal_temp_grid_withdepth_stdev.npy")
combined_temp_grid_withdepth_stdev = np.load("combined_temp_grid_withdepth_stdev.npy")
except:
argo_temp_grid_withdepth_stdev = load_netcdf_grid(with_depth=True, resolution = resolution, mean=argo_temp_grid_withdepth)
seal_temp_grid_withdepth_stdev = load_netcdf_grid("seal", with_depth=True, resolution = resolution, mean=seal_temp_grid_withdepth)
combined_temp_grid_withdepth_stdev = load_netcdf_grid("*", with_depth=True, resolution = resolution, mean=combined_temp_grid_withdepth)
np.save("argo_temp_grid_withdepth_stdev", argo_temp_grid_withdepth_stdev)
np.save("seal_temp_grid_withdepth_stdev", seal_temp_grid_withdepth_stdev)
np.save("combined_temp_grid_withdepth_stdev", combined_temp_grid_withdepth_stdev)
print(argo_temp_grid_withdepth_stdev.shape)
#print(argo_temp_grid_withdepth_stdev[0,0,:],argo_temp_grid_withdepth[0,0,:])
argomean, stdev = [argo_temp_grid_withdepth_stdev[0,0,:], argo_temp_grid_withdepth[0,0,:]]
print(argomean.shape, stdev.shape)
for s in range (0,len(argomean)):
print(s,argomean[s],stdev[s])
(15, 360, 200) (200,) (200,) 0 0.885078899289779 3.406360979433413 1 0.9988122771173085 3.3157297920536353 2 1.0270029350865955 3.386472483476003 3 1.010605817628795 3.3632362650500403 4 1.0252340411131378 3.293358439829812 5 1.01473641644723 3.3054000600179037 6 0.9887591306216763 3.234228641646249 7 0.9679549067768887 3.138942355694978 8 0.9687799560026398 3.066155158298116 9 0.964259996508486 3.046098829994739 10 0.8530076892156978 2.77432834924157 11 0.8186467895643343 2.6918114752009297 12 0.7680624896042024 2.5548307418823244 13 0.7567045228888772 2.540033813250267 14 0.6979674855090975 2.4469851588501648 15 0.6610973533110266 2.388016838138386 16 0.6558932557617471 2.361967638615639 17 0.6896429894718946 2.377774130913519 18 0.6405453780834807 2.330421844497323 19 0.6580049671213198 2.366783275206884 20 0.6098801509317064 2.357095058002169 21 0.6395568438562833 2.3234067084425587 22 0.6448821426620494 2.3149507905616136 23 0.5801722015839754 2.356969274007357 24 0.53267841869067 2.3658066180444535 25 0.4787793549678399 2.473793046227817 26 0.43541491356283796 2.4819062165915966 27 0.37142466738840174 2.5517627505932823 28 0.34738782962804515 2.607932171579135 29 0.339777249736777 2.6501175971592175 30 0.33398208552156294 2.6552541983329645 31 0.30451292267642066 2.689133310317993 32 0.30696012824134306 2.6692257927310084 33 0.3006152472591102 2.6839015835621316 34 0.2752850083694053 2.6960002064704893 35 0.26350202760045743 2.7009692338796762 36 0.2792438224229195 2.689875032220568 37 0.2764682925138239 2.6821965624074466 38 0.26687715472710577 2.6866031260717484 39 0.27010652638444804 2.677758652588417 40 0.25960608725630235 2.662781283259392 41 0.22796535734312656 2.6805423679998364 42 0.1954772437346598 2.6708275408580384 43 0.17187063842758643 2.6710678722898837 44 0.16660867761459625 2.6470793618096247 45 0.16404311269361638 2.6438334306081135 46 0.14723450973243596 2.6379656133980585 47 0.13840563883949064 2.6394753221605645 48 0.14019797995912658 2.6461166103680926 49 0.13646882429368734 2.6629150721986417 50 0.11941947695270799 2.66408345301946 51 0.11950870110326002 2.6627287622225486 52 0.10342887318332163 2.651105148750439 53 0.11611593659816372 2.6524561120752703 54 0.10211777221488016 2.6643548280962053 55 0.0962642230916455 2.66907026474936 56 0.1014485002447104 2.671357044151851 57 0.1004005992401329 2.65904989639918 58 0.09451790242809628 2.6442631294852807 59 0.0832121049512599 2.6218064062057005 60 0.07943381879518607 2.603448275862069 61 0.07553759580545918 2.6021785140037537 62 0.07191198905199282 2.596999950576247 63 0.07612483377884778 2.594158827312409 64 0.06788356469240962 2.5883036894457683 65 0.07465941884064772 2.5788002689679463 66 0.07121779685118802 2.5745423325037553 67 0.07264016221945983 2.570578976681358 68 0.06023895780857734 2.5582679169518605 69 0.061882431153989136 2.5541584453885515 70 0.06011015091764096 2.5438248483758223 71 0.05553506357819721 2.5392361467534847 72 0.05575312301536279 2.5356033876024444 73 0.05652110900742914 2.5273818492889406 74 0.0543692617432526 2.5221378556613265 75 0.06095933960122688 2.516084630610579 76 0.05497390721328046 2.5127855198723927 77 0.05690282765476441 2.513473527473316 78 0.056391123109074334 2.5082406997680664 79 0.05538512439374496 2.5070690204357278 80 0.05734824347638948 2.4909474807873107 81 0.05276936706932022 2.486627077652236 82 0.051102824874305466 2.477232371057783 83 0.049772861851212664 2.4701550746786185 84 0.05444481909010995 2.462055656645033 85 0.05311064944090903 2.4514736585449755 86 0.050032877638898535 2.446125166756766 87 0.04915043920937369 2.4406331181526184 88 0.05086844597521709 2.4333519052576134 89 0.050781428569792314 2.4291401829635886 90 0.05217033046108383 2.421300152937571 91 0.05064196210666983 2.4168750814029147 92 0.04937789767724743 2.408833450741238 93 0.049160670772637055 2.4072835167249043 94 0.04944694395948274 2.3986665584422924 95 0.05156224270902831 2.3956066311382855 96 0.05101647434187872 2.3914184657010167 97 0.0512445469188101 2.385647981255143 98 0.04715627175230448 2.3805000952311923 99 0.04604398315480987 2.3764262003976793 100 0.050301650193789384 2.3670878326683713 101 0.04632860500959619 2.3638948139391447 102 0.046527801057915846 2.3570000778545035 103 0.04673741890937641 2.348611072257713 104 0.04640289753485496 2.3456948692515747 105 0.04686031576517745 2.3344034688514577 106 0.04332060409382052 2.330553582736424 107 0.04408705540510528 2.323309222134677 108 0.04535004553506316 2.320625058242253 109 0.044078446428925444 2.315271296743619 110 0.04377884899000489 2.3044136968152276 111 0.0407608075378575 2.299357150282179 112 0.03930798304617833 2.294070093255294 113 0.04116365575587148 2.2878703188013145 114 0.04394590097210776 2.2871929762656227 115 0.042232092979545895 2.277175229892396 116 0.0396330717462302 2.2739123461539283 117 0.03975091362139861 2.267839193344116 118 0.039643815466092705 2.2632142645972118 119 0.043401896148458075 2.256524547201688 120 0.042123197742706335 2.24931471436112 121 0.03973026871772996 2.242600024830211 122 0.03956290229640913 2.238779585240251 123 0.04091210421523333 2.2317038068064936 124 0.041357055721067154 2.2287929386928163 125 0.04130288648577337 2.2173272609710692 126 0.04004833537462321 2.216912369979055 127 0.03871143896160243 2.208527183532715 128 0.039499057709008016 2.2024814641034163 129 0.04082841409698702 2.1993226582004177 130 0.04314982166523841 2.1888571637017384 131 0.041543104926593286 2.1855087782207288 132 0.044689432703738936 2.178618218682029 133 0.043843147074903686 2.171818178350275 134 0.04574259371825896 2.1665254204960194 135 0.04325069252954223 2.1602240102044465 136 0.043070572125732266 2.152678574834551 137 0.04409696653033099 2.1472223069932728 138 0.04577291853362123 2.143543900104991 139 0.04514203950625671 2.137947287475854 140 0.044560991441594075 2.128338890560603 141 0.044461252019087585 2.122267872095108 142 0.04516149473163497 2.1162857839039395 143 0.04660426837701246 2.1089259253607855 144 0.04922777584676679 2.1066205748196305 145 0.04833361751216056 2.0947017627849913 146 0.047354701744565054 2.090309056368741 147 0.04710193817948876 2.0840703018924644 148 0.04918936509126414 2.0774546016346322 149 0.049569630746547653 2.0740862098233452 150 0.05042300966849336 2.0645667652289075 151 0.04929160892068561 2.057641533185851 152 0.049014535265846426 2.052418171275746 153 0.04742252165398071 2.0461355872073415 154 0.04966537351281795 2.0406206599597274 155 0.051124368037648076 2.0300174470533405 156 0.05019612393251732 2.0288771913762678 157 0.05031346424648784 2.0200944036807655 158 0.04941821016882442 2.0121999567205258 159 0.052521692345925755 2.0063869933928213 160 0.0495863892924325 1.9978928651128496 161 0.049979994715309034 1.9920000249689276 162 0.052255554981481565 1.9862320678574699 163 0.053193624241901744 1.9766481187608507 164 0.05329419595241589 1.973389791230024 165 0.05210199190902407 1.9617903924757434 166 0.054313242818351926 1.9527307198597834 167 0.0546730890965354 1.9471785851887293 168 0.05563153616300946 1.9389272299679843 169 0.054955892749498754 1.9337543500097174 170 0.054881703885290734 1.92312071652248 171 0.051110494579914886 1.9175333062807718 172 0.05478543225403508 1.9056923572833722 173 0.05327904546494806 1.9009465319769723 174 0.05614157806685546 1.8973508634065326 175 0.054814815708991466 1.8856609514204121 176 0.055368131626691554 1.8784544684670188 177 0.05446607267921217 1.8717759559894431 178 0.056872926887882126 1.8626295548898202 179 0.059053950612033765 1.8577719445814167 180 0.060619991327868994 1.8479833245277404 181 0.058872956460463206 1.838222274073848 182 0.05894663167092525 1.834285821233477 183 0.0580443527369214 1.8272713018675981 184 0.06251718825952234 1.8202364054593172 185 0.062248188417342636 1.8080862123390724 186 0.062434452040562465 1.8067858197859354 187 0.0599354587003413 1.7953332971643519 188 0.06078765928460939 1.7885000970628526 189 0.058005064780800926 1.7847376768706276 190 0.06210909088302475 1.7716363668441772 191 0.05922835907187037 1.7675357035228185 192 0.05928831301635136 1.7579055642181973 193 0.0579684071145023 1.7500536101205009 194 0.06403721604853121 1.738946482539177 195 0.057441273485671865 1.7328102896953452 196 0.059263016386028634 1.722685213442202 197 0.059560802064034024 1.7129635550759055 198 0.05356288142892028 1.712132157019849 199 0.044638491914242316 1.7177334785461427
# Load / Create the Lat-Long-Depth Grid of the Climatological Salinity Field
resolution = 1
try:
argo_sal_grid_withdepth = np.load("argo_sal_grid_withdepth.npy")
seal_sal_grid_withdepth = np.load("seal_sal_grid_withdepth.npy")
combined_sal_grid_withdepth = np.load("combined_sal_grid_withdepth.npy")
except:
argo_sal_grid_withdepth = load_netcdf_grid(target_var = "sal", with_depth=True, resolution = resolution)
seal_sal_grid_withdepth = load_netcdf_grid("seal", target_var = "sal", with_depth=True, resolution = resolution)
combined_sal_grid_withdepth = np.nanmean((argo_sal_grid_withdepth, seal_sal_grid_withdepth), axis=0)
np.save("argo_sal_grid_withdepth", argo_sal_grid_withdepth)
np.save("seal_sal_grid_withdepth", seal_sal_grid_withdepth)
np.save("combined_sal_grid_withdepth", combined_sal_grid_withdepth)
print(argo_temp_grid_withdepth.shape, seal_temp_grid_withdepth.shape, combined_temp_grid_withdepth.shape)
(15, 360, 200) (15, 360, 200) (15, 360, 200)
#dbars = [0, 4, 9, 14, 19, 29, 49, 69, 99, 124, 149, 174, 189, 199]
dbars = [0]
# Plot Climatological Temperature Field
tmin, tmax = [np.nanmin(combined_temp_grid_withdepth[:,:,db]),np.nanmax(combined_temp_grid_withdepth[:,:,db])]
tmin, tmax = [np.rint(tmin), np.ceil(tmax)]
for db in dbars:
#Argo data\n0-10 decibar mean temperature for {}\nMean across all points = {}".format(calendar.month_abbr[month], argomean)
#dbs = "{}-{} dbar".format(db*10, (db+1)*10)
dbs = "between {}-{} dbar".format(db*10, (db+1)*10)
#plot_grid(argo_temp_grid_withdepth[:,:,db], vmin=tmin, vmax=tmax)
#plot_grid(seal_temp_grid_withdepth[:,:,db], vmin=tmin, vmax=tmax)
print("valor minimo e maximo", np.nanmin(combined_temp_grid_withdepth[:,:,db]),np.nanmax(combined_temp_grid_withdepth[:,:,db]))
plot_grid(combined_temp_grid_withdepth[:,:,db], resolution=resolution, vmin=tmin, vmax=tmax)
#plot_grid(argo_temp_grid_withdepth[:,:,db], title=a_title + "\n" + mtf + dbs, resolution=resolution, vmin=tmin, vmax=tmax)
#plot_grid(seal_temp_grid_withdepth[:,:,db], title=s_title + "\n" + mtf + dbs, resolution=resolution, vmin=tmin, vmax=tmax)
#plot_grid(combined_temp_grid_withdepth[:,:,db], title=comb_title+ "\n"+ mtf + dbs, resolution=resolution, vmin=tmin, vmax=tmax)
#plot_grid(combined_temp_grid_withdepth[:,:,db], title=dbs, resolution=resolution, vmin=tmin, vmax=tmax)
#interp = interp_nans(combined_temp_grid_withdepth[:,:,db])
#plot_grid(interp, title=interp_title + "\n" + mtf + dbs, resolution=resolution, vmin=tmin, vmax=tmax)
valor minimo e maximo -1.9446429184504919 5.1874998807907104
# Plot Climatological Salinity Field
smin, smax = [np.nanmin(combined_sal_grid_withdepth[:,:,db]),np.nanmax(combined_sal_grid_withdepth[:,:,db])]
smin, smax = [np.floor(smin), np.ceil(smax)]
#smin, smax = [33.8,35];
for db in dbars:
#Argo data\n0-10 decibar mean temperature for {}\nMean across all points = {}".format(calendar.month_abbr[month], argomean)
#dbs = "{}-{} dbar".format(db*10, (db+1)*10)
dbs = "between {}-{} dbar".format(db*10, (db+1)*10)
#plot_grid(argo_sal_grid_withdepth[:,:,db], title=a_title + "\n" + msf + dbs , cbtitle="Salinity", resolution=resolution, vmin=smin, vmax=smax)
#plot_grid(seal_sal_grid_withdepth[:,:,db], title=s_title + "\n" + msf + dbs , cbtitle="Salinity", resolution=resolution, vmin=smin, vmax=smax)
#plot_grid(combined_sal_grid_withdepth[:,:,db], title=comb_title + "\n" + msf + dbs , cbtitle="Salinity", resolution=resolution, vmin=smin, vmax=smax)
#interp = interp_nans(combined_sal_grid_withdepth[:,:,db])
#plot_grid(interp, title=interp_title + "\n" + msf + dbs, cbtitle="Salinity", resolution=resolution, vmin=smin, vmax=smax)
print("valor minimo e maximo", np.nanmin(combined_sal_grid_withdepth[:,:,db]),np.nanmax(combined_sal_grid_withdepth[:,:,db]))
plot_grid(combined_sal_grid_withdepth[:,:,db], cbtitle="Salinity", resolution=resolution, vmin=smin, vmax=smax)
valor minimo e maximo 32.604942321777344 34.585634867350265
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
#m.drawcoastlines()
#m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
m.etopo()
plt.title("Bathymetry south of 60S")
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Text(0.5, 1.0, 'Bathymetry south of 60S')
import pylab as plot
params = {'legend.fontsize': 14,
'legend.handlelength': 1,
'legend.labelspacing': 0.25,
'legend.handletextpad': 0.2,
'legend.frameon': False,
'legend.markerscale': 3.0}
plot.rcParams.update(params)
def format_axes_b(fig):
for i, ax in enumerate(fig.axes):
ax.text(-1.5, 1950, "4.%d" % (i+1), va="center", ha="center", fontsize=14, color='black')
ax.tick_params(axis="both", direction='inout', which='major', labelsize=16, labelbottom=True, labelleft=True)
ax.set_ylim(p_max, p_min);
ax.set_xlim(t_min, t_max)
ax1.set_ylabel("Depth (dbar)", fontsize="18")
ax2.set_ylabel("Depth (dbar)", fontsize="18")
ax2.set_xlabel("Temperature", fontsize="18")
ax6.set_xlabel("Temperature", fontsize="18")
ax1.tick_params(labelbottom=False)
i = fig = b = c = d = e = f = g = long = 0
resolution = 1
cbtitle = title = parameter = "Temperature ºC"
def plot_cross_slope_(grid1, grid2, title=title, lon=long, i = i, fig = fig, cbtitle = cbtitle, resolution = resolution,
labelleft=True, labelbottom=True, tmin=b, tmax=c, smin=d, smax=e, dmin = f, dmax = g):
ax = fig.add_subplot(2, 2, i)
ax.tick_params(axis="both", direction='inout', labelsize=14, labelbottom=True, labelleft=True, labeltop=False)
y = np.arange(-60, -75, -1/resolution)
z = np.arange(0, 2000, 10)
levels = np.arange(np.nanmin(grid1), np.nanmax(grid1), .01)
im = ax.contourf(y, z, grid1.T, cmap="nipy_spectral", levels = levels, extend="both")
ax.text(-60.5, 1900, "(8.%d) {}º".format(lon) % (i), va="center", ha="left", fontsize=12, color='black')
#ax.set_title("{} at {} W at {} resolution".format(title,lon, 1 / resolution))
#ax.set_title("at {}".format(lon))
cbticks = np.arange(-2.5,tmax,((tmax-tmin)/9))
cb = fig.colorbar(im, pad=0.025, format='%.2f', extend='both',ticks=cbticks)
cb_ticklabel_size = 14 # Adjust as appropriate.
cb.ax.tick_params(labelsize=cb_ticklabel_size)
label_size = 16 # Adjust as appropriate.
#ax.tick_params(axis="both", direction='inout', which='major')
if i<3:
ax.tick_params(labelbottom=False)
ax.tick_params(labeltop=False, direction='inout')
else:
ax.tick_params(labelbottom=True)
ax.set_xlabel("° Latitude", fontsize=label_size)
if (i % 2) == 0:
ax.tick_params(labelleft=False)
cb.set_label(cbtitle, rotation=270, fontsize="16", labelpad=12.5)
else:
ax.set_ylabel("Depth (dbar)", fontsize=label_size)
smin, smax = [33, 35]
contour_levels = [33.3, 33.5, 33.7, 33.9, 34.1, 34.3, 34.5, 34.7, 34.8]
#contour_levels = np.linspace(smin, smax, num=21) # Adjust as appropriate.
print(contour_levels)
CS = ax.contour(y, z, grid2.T, #10,
colors='k', vmin=smin, vmax=smax, levels=contour_levels)
ax.clabel(CS, fontsize=9, inline=1, fmt='%1.2f')
#cb = fig.colorbar(im, CS, pad=0.025, format='%.2f', extend='both',ticks=cbticks)
db = oceansdb.ETOPO()
y = np.arange(-60, -75, -.01)
h = -db['topography'].extract(lat=y, lon=long)["height"]
#ax.plot(y, h)
ax.fill_between(y, 5000, h, color='black')
ax.set_ylim(2000, 0)
ax.set_xlim(-60, -74.5)
#longit = (0, - 178.5, -90.5, 0.5, 90) # Half Degree Resolution
longit = (0,180,-90,0.,90) # One Degree Resolution
print(longit[:])
z_grid = (0, 0, 90, 180, 270) # One Degree Resolution
#z_grid = (0, 3, 179, 361, 540) # Half Degree Resolution
print(z_grid[:])
fig = plt.figure(figsize=(25,10))
fig.subplots_adjust(hspace=0.05, wspace=0.005)
for i in range(1, 5):
a = z_grid[i]
long = longit[i]
print(a,long)
xsection_temp = combined_temp_grid_withdepth[:,a,:]
xsection_temp = interp_nans(xsection_temp)
xsection_sal = argo_sal_grid_withdepth[:,a,:]
xsection_sal = interp_nans(xsection_sal)
b,c,d,e = [np.nanmin(xsection_temp),np.nanmax(xsection_temp),np.nanmin(xsection_sal),np.nanmax(xsection_sal)]
print(b,c,d,e)
plot_cross_slope_(xsection_temp, xsection_sal, lon = long, i = i, fig = fig, tmin=b, tmax=c, smin=d, smax=e)
#fig.suptitle("Cross section of {} at {} resolution".format(title, 1 / resolution), fontsize="20")
plt.show()
(0, 180, -90, 0.0, 90) (0, 0, 90, 180, 270) 0 180 -1.7390999794006348 3.3632362650500403 33.731198120117185 34.73633130391439 [33.3, 33.5, 33.7, 33.9, 34.1, 34.3, 34.5, 34.7, 34.8] 90 -90 -1.8242663070559502 3.86983315149943 33.45185068491343 34.73500061035156 [33.3, 33.5, 33.7, 33.9, 34.1, 34.3, 34.5, 34.7, 34.8] 180 0.0 -1.7166638456541916 1.002409734641369 33.870741965553975 34.69793701171875 [33.3, 33.5, 33.7, 33.9, 34.1, 34.3, 34.5, 34.7, 34.8] 270 90 -1.8788625001907349 1.8338245834623064 33.369389257123395 34.735844135284424 [33.3, 33.5, 33.7, 33.9, 34.1, 34.3, 34.5, 34.7, 34.8]
longitude = long = 180
cbtitle = title = "Temperature (°C)"
def plot_cross_slope(grid, title=title, lon=long, cbtitle=title, resolution = resolution):
fig, ax = plt.subplots(1, 1, figsize=(10,5))
y = np.arange(-60, -75, -1/resolution)
z = np.arange(0, 2000, 10)
levels = np.arange(np.nanmin(grid), np.nanmax(grid), .01)
im = ax.contourf(y, z, grid.T, cmap="nipy_spectral", levels = levels, extend="both")
ax.set_xlabel("° Latitude")
ax.set_ylabel("Depth (dbar)")
ax.set_title("{} at {} W at {} resolution".format(title,longitude, 1 / resolution))
cb = fig.colorbar(im, pad=0.025)
cb.set_label(cbtitle, rotation=270, fontsize="14", labelpad=12.5)
#minor_thresholds=(2, 0.5)
db = oceansdb.ETOPO()
y = np.arange(-60, -75, -.01)
h = -db['topography'].extract(lat=y, lon=long)["height"]
#ax.plot(y, h)
ax.fill_between(y, 5000, h, color='black')
ax.set_ylim(2000, 0)
ax.set_xlim(-60, -74.5)
plt.show()
zonal_grid = 90
filtered = combined_temp_grid_withdepth[:,zonal_grid,:]
interp = interp_nans(filtered)
plot_cross_slope(filtered, resolution = resolution)
plot_cross_slope(interp, resolution = resolution)
def dens_smow(T):
'''
Function calculates density of standard mean ocean water (pure water) using EOS 1980.
INPUT: T = temperature [degree C (ITS-90)]
OUTPUT: dens = density [kg/m^3]
'''
a0 = 999.842594;a1 = 6.793952e-2;a2 = -9.095290e-3;a3 = 1.001685e-4;a4 = -1.120083e-6;a5 = 6.536332e-9
T68 = T * 1.00024
dens = (a0 + (a1 + (a2 + (a3 + (a4 + a5*T68)*T68)*T68)*T68)*T68)
return dens
def dens0(S, T):
'''
Function calculates density of sea water at atmospheric pressure
USAGE: dens0 = dens0(S,T)
DESCRIPTION:
Density of Sea Water at atmospheric pressure using
UNESCO 1983 (EOS 1980) polynomial.
INPUT: (all must have same dimensions)
S = salinity [psu (PSS-78)]
T = temperature [degree C (ITS-90)]
OUTPUT:
dens0 = density [kg/m^3] of salt water with properties S,T,
P=0 (0 db gauge pressure)
'''
assert S.shape == T.shape
T68 = T * 1.00024
# UNESCO 1983 eqn(13) p17.
b0 = 8.24493e-1;b1 = -4.0899e-3;b2 = 7.6438e-5;b3 = -8.2467e-7;b4 = 5.3875e-9
c0 = -5.72466e-3;c1 = +1.0227e-4;c2 = -1.6546e-6
d0 = 4.8314e-4
dens = (dens_smow(T) + (b0+ (b1+(b2+(b3+b4*T68)*T68)*T68)*T68)*S + (c0+(c1+c2*T68)*T68)*S*np.sqrt(S)+ d0*(S**2) )
return dens
parameter = "Density (kg/m3)"
dens = dens0(combined_sal_grid_withdepth, combined_temp_grid_withdepth)
sigma = dens[:,:,:] - 1000
print(dens.shape, sigma.shape)
for s in range (0,len(dens)):
print(s,dens[s],sigma[s])
(15, 360, 200) (15, 360, 200)
0 [[1026.99899305 1027.00295946 1027.00104598 ... 1027.78233243
1027.78232547 1027.78183647]
[1027.0299895 1027.02664621 1027.02591856 ... 1027.78432079
1027.78499903 1027.78467147]
[1027.06583993 1027.03544742 1027.03972517 ... 1027.78706838
1027.78761044 1027.78811721]
...
[1027.0231622 1027.02920716 1027.02812946 ... 1027.78380222
1027.78386284 1027.78165242]
[1026.99677127 1026.99427819 1026.99312125 ... 1027.78529784
1027.78531989 1027.7822855 ]
[1027.06763691 1027.05785841 1027.05293542 ... 1027.77933789
1027.78013152 1027.7809418 ]] [[26.99899305 27.00295946 27.00104598 ... 27.78233243 27.78232547
27.78183647]
[27.0299895 27.02664621 27.02591856 ... 27.78432079 27.78499903
27.78467147]
[27.06583993 27.03544742 27.03972517 ... 27.78706838 27.78761044
27.78811721]
...
[27.0231622 27.02920716 27.02812946 ... 27.78380222 27.78386284
27.78165242]
[26.99677127 26.99427819 26.99312125 ... 27.78529784 27.78531989
27.7822855 ]
[27.06763691 27.05785841 27.05293542 ... 27.77933789 27.78013152
27.7809418 ]]
1 [[1027.08762482 1027.06330596 1027.07063059 ... 1027.78140186
1027.78359607 1027.78588772]
[1027.10533451 1027.09862571 1027.09705175 ... 1027.79376575
1027.7942402 1027.79142985]
[1027.1032459 1027.08680705 1027.08828786 ... 1027.79619178
1027.79705787 1027.79659641]
...
[1026.98646 1026.98176037 1026.98437123 ... 1027.78473021
1027.78579028 1027.78607999]
[1027.02329228 1027.02236432 1027.02517175 ... 1027.78684211
1027.78705131 1027.78501985]
[1026.99358329 1026.99960295 1027.00135489 ... 1027.78801106
1027.78699075 1027.78565308]] [[27.08762482 27.06330596 27.07063059 ... 27.78140186 27.78359607
27.78588772]
[27.10533451 27.09862571 27.09705175 ... 27.79376575 27.7942402
27.79142985]
[27.1032459 27.08680705 27.08828786 ... 27.79619178 27.79705787
27.79659641]
...
[26.98646 26.98176037 26.98437123 ... 27.78473021 27.78579028
27.78607999]
[27.02329228 27.02236432 27.02517175 ... 27.78684211 27.78705131
27.78501985]
[26.99358329 26.99960295 27.00135489 ... 27.78801106 27.78699075
27.78565308]]
2 [[1027.05744738 1027.04915021 1027.04853385 ... 1027.79552049
1027.79563162 1027.79537031]
[1027.0756876 1027.06289491 1027.06251589 ... 1027.7972256
1027.79779489 1027.79806875]
[1027.07759048 1027.05093199 1027.05236452 ... 1027.79955374
1027.80019021 1027.79881744]
...
[1026.99787548 1026.98903185 1026.9875189 ... 1027.79749658
1027.79811166 1027.79932697]
[1027.01037535 1026.99789388 1027.00243965 ... 1027.79875745
1027.79793765 1027.79825587]
[1027.01105427 1027.00157867 1027.00347931 ... 1027.79752435
1027.7973029 1027.79733141]] [[27.05744738 27.04915021 27.04853385 ... 27.79552049 27.79563162
27.79537031]
[27.0756876 27.06289491 27.06251589 ... 27.7972256 27.79779489
27.79806875]
[27.07759048 27.05093199 27.05236452 ... 27.79955374 27.80019021
27.79881744]
...
[26.99787548 26.98903185 26.9875189 ... 27.79749658 27.79811166
27.79932697]
[27.01037535 26.99789388 27.00243965 ... 27.79875745 27.79793765
27.79825587]
[27.01105427 27.00157867 27.00347931 ... 27.79752435 27.7973029
27.79733141]]
3 [[1027.04547284 1027.04852503 1027.05116031 ... 1027.80477258
1027.80492471 1027.80444118]
[1027.06618294 1027.06629816 1027.06792448 ... 1027.80535284
1027.80473182 1027.80579807]
[1027.00735897 1027.01490862 1027.01932091 ... 1027.8014453
1027.80211524 1027.80154763]
...
[1026.98340598 1026.97932146 1026.98088378 ... 1027.80523537
1027.80682714 1027.80467324]
[1026.97879385 1026.98526469 1026.98915093 ... 1027.80513169
1027.8053845 1027.80576159]
[1026.97739416 1026.98874109 1026.99219919 ... 1027.80432304
1027.80480517 1027.80397517]] [[27.04547284 27.04852503 27.05116031 ... 27.80477258 27.80492471
27.80444118]
[27.06618294 27.06629816 27.06792448 ... 27.80535284 27.80473182
27.80579807]
[27.00735897 27.01490862 27.01932091 ... 27.8014453 27.80211524
27.80154763]
...
[26.98340598 26.97932146 26.98088378 ... 27.80523537 27.80682714
27.80467324]
[26.97879385 26.98526469 26.98915093 ... 27.80513169 27.8053845
27.80576159]
[26.97739416 26.98874109 26.99219919 ... 27.80432304 27.80480517
27.80397517]]
4 [[1027.27705105 1027.27810375 1027.28326856 ... 1027.81720307
1027.81746304 1027.82057178]
[1027.18538306 1027.19881172 1027.21658779 ... 1027.8133228
1027.8151226 1027.81406861]
[1027.29032961 1027.30077378 1027.31484944 ... 1027.8133645
1027.81581038 1027.81436341]
...
[1027.10852631 1027.12314614 1027.13407602 ... 1027.81880218
1027.81834026 1027.81957225]
[1027.17943537 1027.17565245 1027.18391111 ... 1027.8191804
1027.81968879 1027.81925067]
[1027.20358783 1027.21158268 1027.21962103 ... 1027.81801609
1027.8187011 1027.81999676]] [[27.27705105 27.27810375 27.28326856 ... 27.81720307 27.81746304
27.82057178]
[27.18538306 27.19881172 27.21658779 ... 27.8133228 27.8151226
27.81406861]
[27.29032961 27.30077378 27.31484944 ... 27.8133645 27.81581038
27.81436341]
...
[27.10852631 27.12314614 27.13407602 ... 27.81880218 27.81834026
27.81957225]
[27.17943537 27.17565245 27.18391111 ... 27.8191804 27.81968879
27.81925067]
[27.20358783 27.21158268 27.21962103 ... 27.81801609 27.8187011
27.81999676]]
5 [[1027.2999137 1027.32267489 1027.33265384 ... 1027.82294859
1027.82363688 1027.82455419]
[1027.23794141 1027.2659819 1027.28145905 ... 1027.82099785
1027.82124859 1027.82094525]
[1027.28474753 1027.29771077 1027.31396759 ... 1027.82038928
1027.82189292 1027.82337535]
...
[1027.16918847 1026.98965465 1027.07551569 ... 1027.82737721
1027.82764896 1027.82788357]
[1027.09855568 1027.00668253 1027.11093043 ... 1027.82186317
1027.82194679 1027.82334531]
[1027.06562493 1027.08865289 1027.09559232 ... 1027.82159983
1027.82131516 1027.8220611 ]] [[27.2999137 27.32267489 27.33265384 ... 27.82294859 27.82363688
27.82455419]
[27.23794141 27.2659819 27.28145905 ... 27.82099785 27.82124859
27.82094525]
[27.28474753 27.29771077 27.31396759 ... 27.82038928 27.82189292
27.82337535]
...
[27.16918847 26.98965465 27.07551569 ... 27.82737721 27.82764896
27.82788357]
[27.09855568 27.00668253 27.11093043 ... 27.82186317 27.82194679
27.82334531]
[27.06562493 27.08865289 27.09559232 ... 27.82159983 27.82131516
27.8220611 ]]
6 [[1027.1384504 1027.12491052 1027.09047823 ... 1027.82869667
1027.82940302 1027.83024616]
[1027.1160797 1027.12850828 1027.14063806 ... 1027.82449429
1027.82422621 1027.82470347]
[1027.18978392 1027.21088337 1027.22154819 ... 1027.82414413
1027.82472323 1027.82727087]
...
[1027.24347343 1027.38880026 1027.37080296 ... 1027.83625994
1027.83641888 1027.83717489]
[1027.27737227 1027.37700948 1027.34697339 ... 1027.83693295
1027.83706552 1027.83659592]
[1027.36937393 1027.39320841 1027.39361969 ... 1027.8347176
1027.83496188 1027.83402764]] [[27.1384504 27.12491052 27.09047823 ... 27.82869667 27.82940302
27.83024616]
[27.1160797 27.12850828 27.14063806 ... 27.82449429 27.82422621
27.82470347]
[27.18978392 27.21088337 27.22154819 ... 27.82414413 27.82472323
27.82727087]
...
[27.24347343 27.38880026 27.37080296 ... 27.83625994 27.83641888
27.83717489]
[27.27737227 27.37700948 27.34697339 ... 27.83693295 27.83706552
27.83659592]
[27.36937393 27.39320841 27.39361969 ... 27.8347176 27.83496188
27.83402764]]
7 [[1027.26458045 1027.4203274 1027.27954627 ... 1027.83574389
1027.83603589 1027.83758727]
[1027.26012975 1027.31461797 1027.31477819 ... 1027.83177405
1027.83352136 1027.82913869]
[1027.23647799 1027.27899185 1027.27036331 ... 1027.83439303
1027.83442012 1027.83428168]
...
[1027.30292671 1027.40533612 1027.42166955 ... 1027.83447207
1027.8345736 1027.83445903]
[1027.33296155 1027.391663 1027.38844008 ... 1027.83435194
1027.83447538 1027.83457003]
[1027.26808569 1027.36350799 1027.2869419 ... 1027.83444138
1027.83430328 nan]] [[27.26458045 27.4203274 27.27954627 ... 27.83574389 27.83603589
27.83758727]
[27.26012975 27.31461797 27.31477819 ... 27.83177405 27.83352136
27.82913869]
[27.23647799 27.27899185 27.27036331 ... 27.83439303 27.83442012
27.83428168]
...
[27.30292671 27.40533612 27.42166955 ... 27.83447207 27.8345736
27.83445903]
[27.33296155 27.391663 27.38844008 ... 27.83435194 27.83447538
27.83457003]
[27.26808569 27.36350799 27.2869419 ... 27.83444138 27.83430328
nan]]
8 [[1027.05506344 1027.06897555 1027.03980098 ... nan
nan nan]
[1027.2243647 1027.23316346 1027.22958762 ... 1027.8280465
nan nan]
[1027.22132071 1027.23804441 1027.23680329 ... nan
nan nan]
...
[1027.10657307 1027.40729612 1027.43319049 ... nan
nan 1027.82930458]
[1026.96728893 nan 1026.95306711 ... nan
nan nan]
[1027.02632216 1027.00176167 1027.02197893 ... nan
nan nan]] [[27.05506344 27.06897555 27.03980098 ... nan nan
nan]
[27.2243647 27.23316346 27.22958762 ... 27.8280465 nan
nan]
[27.22132071 27.23804441 27.23680329 ... nan nan
nan]
...
[27.10657307 27.40729612 27.43319049 ... nan nan
27.82930458]
[26.96728893 nan 26.95306711 ... nan nan
nan]
[27.02632216 27.00176167 27.02197893 ... nan nan
nan]]
9 [[1027.08999673 1027.11083709 1027.09727198 ... nan
nan 1027.83211669]
[1027.13613265 1027.11245518 1027.12425021 ... nan
nan nan]
[1027.14749407 1027.15779564 1027.23757442 ... 1027.82708727
nan 1027.82733602]
...
[1027.15974221 1027.32634585 1027.30202916 ... nan
nan 1027.83308676]
[1027.33333924 1027.34744194 1027.34250183 ... nan
nan 1027.83270953]
[1027.02805687 1027.08605662 1027.10581737 ... nan
nan 1027.82894134]] [[27.08999673 27.11083709 27.09727198 ... nan nan
27.83211669]
[27.13613265 27.11245518 27.12425021 ... nan nan
nan]
[27.14749407 27.15779564 27.23757442 ... 27.82708727 nan
27.82733602]
...
[27.15974221 27.32634585 27.30202916 ... nan nan
27.83308676]
[27.33333924 27.34744194 27.34250183 ... nan nan
27.83270953]
[27.02805687 27.08605662 27.10581737 ... nan nan
27.82894134]]
10 [[1027.13182055 1027.12600232 1027.14167341 ... 1027.82949078
nan 1027.83170932]
[1027.00265202 1027.38426271 1027.07909986 ... 1027.82997645
nan 1027.83223476]
[1026.86206623 1027.01394439 1027.05069811 ... 1027.82952679
nan 1027.83125535]
...
[ nan nan nan ... nan
nan nan]
[ nan nan nan ... nan
nan nan]
[1026.98676898 1026.98682183 1026.98824205 ... nan
nan 1027.83138017]] [[27.13182055 27.12600232 27.14167341 ... 27.82949078 nan
27.83170932]
[27.00265202 27.38426271 27.07909986 ... 27.82997645 nan
27.83223476]
[26.86206623 27.01394439 27.05069811 ... 27.82952679 nan
27.83125535]
...
[ nan nan nan ... nan nan
nan]
[ nan nan nan ... nan nan
nan]
[26.98676898 26.98682183 26.98824205 ... nan nan
27.83138017]]
11 [[ nan nan nan ... nan
nan nan]
[1026.78843318 nan nan ... nan
nan nan]
[1027.20248915 1027.19803927 1027.20523781 ... nan
nan 1027.83360126]
...
[1027.30936367 1027.40168722 1027.3474926 ... nan
nan 1027.83619299]
[1027.20018566 1027.56915734 1027.26056306 ... nan
nan 1027.83772627]
[1027.15614646 1027.2355082 1027.25187562 ... nan
nan 1027.83940837]] [[ nan nan nan ... nan nan
nan]
[26.78843318 nan nan ... nan nan
nan]
[27.20248915 27.19803927 27.20523781 ... nan nan
27.83360126]
...
[27.30936367 27.40168722 27.3474926 ... nan nan
27.83619299]
[27.20018566 27.56915734 27.26056306 ... nan nan
27.83772627]
[27.15614646 27.2355082 27.25187562 ... nan nan
27.83940837]]
12 [[1026.91942121 1026.91363808 1026.91384344 ... nan
nan nan]
[1026.83578115 nan nan ... nan
nan nan]
[1027.08688458 1027.37547912 1027.11862432 ... nan
nan nan]
...
[1027.28070589 1027.43362492 1027.43370845 ... nan
nan nan]
[1027.1576939 1027.15990378 1027.16342706 ... nan
nan nan]
[1027.31762706 1027.31875436 1027.32050198 ... nan
nan 1027.83879835]] [[26.91942121 26.91363808 26.91384344 ... nan nan
nan]
[26.83578115 nan nan ... nan nan
nan]
[27.08688458 27.37547912 27.11862432 ... nan nan
nan]
...
[27.28070589 27.43362492 27.43370845 ... nan nan
nan]
[27.1576939 27.15990378 27.16342706 ... nan nan
nan]
[27.31762706 27.31875436 27.32050198 ... nan nan
27.83879835]]
13 [[1027.47882743 1027.59301782 1027.59596959 ... nan
nan 1027.84156357]
[1027.47141839 1027.57543284 1027.57609519 ... nan
nan 1027.83936002]
[1027.45159529 1027.48107557 1027.49083473 ... nan
nan nan]
...
[1027.28147521 1027.2828009 1027.28359764 ... nan
nan nan]
[1027.40087677 1027.40311969 1027.40328682 ... nan
nan nan]
[1027.2921766 1027.30005393 1027.30800155 ... nan
nan 1027.84055086]] [[27.47882743 27.59301782 27.59596959 ... nan nan
27.84156357]
[27.47141839 27.57543284 27.57609519 ... nan nan
27.83936002]
[27.45159529 27.48107557 27.49083473 ... nan nan
nan]
...
[27.28147521 27.2828009 27.28359764 ... nan nan
nan]
[27.40087677 27.40311969 27.40328682 ... nan nan
nan]
[27.2921766 27.30005393 27.30800155 ... nan nan
27.84055086]]
14 [[1027.51420019 1027.54740894 1027.58447019 ... nan
nan nan]
[1027.45304472 1027.53075846 1027.60196157 ... nan
nan nan]
[ nan 1027.64797799 1027.64980316 ... nan
nan nan]
...
[1027.55013912 1027.56937547 1027.61610347 ... nan
nan nan]
[1027.47636282 1027.47916789 1027.47971514 ... nan
nan nan]
[1027.7158472 1027.71784705 1027.72309246 ... nan
nan nan]] [[27.51420019 27.54740894 27.58447019 ... nan nan
nan]
[27.45304472 27.53075846 27.60196157 ... nan nan
nan]
[ nan 27.64797799 27.64980316 ... nan nan
nan]
...
[27.55013912 27.56937547 27.61610347 ... nan nan
nan]
[27.47636282 27.47916789 27.47971514 ... nan nan
nan]
[27.7158472 27.71784705 27.72309246 ... nan nan
nan]]
# SIGMA ==> Density - 1000
parameter = "Sigma \u03C3"
dens = dens0(combined_sal_grid_withdepth, combined_temp_grid_withdepth)
sigma = dens[:,:,:] - 1000
dbars = [50]
for db in dbars:
rhomin, rhomax = [np.nanmin(sigma[:,:,db]),np.nanmax(sigma[:,:,db])]
dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
plot_grid(sigma[:,:,db], title=comb_title + "\n" + msd + dbs,
cbtitle=parameter, resolution=resolution, vmin=rhomin, vmax=rhomax)
#interp = interp_nans(dens[:,:,db])
#print(np.nanmin(interp), np.nanmax(interp))
#plot_grid(interp, title="Interpolated gridded seal+argo density data " + dbs, cbtitle="Density", resolution=resolution, vmin=rhomin, vmax=rhomax)
i = fig = b = c = d = e = f = g = long = 0
resolution = 1
cbtitle = title = parameter = "Density (kg/m3)"
def plot_cross_slope_dens(grid1, grid2, title=title, lon=long, i = i, fig = fig, cbtitle = cbtitle, resolution = resolution,
labelleft=True, labelbottom=True, tmin=b, tmax=c, smin=d, smax=e, dmin = f, dmax = g):
ax = fig.add_subplot(2, 2, i)
ax.tick_params(axis="both", direction='inout', labelsize=14, labelbottom=True, labelleft=True, labeltop=False)
y = np.arange(-60, -75, -1/resolution)
z = np.arange(0, 2000, 10)
levels = np.arange(np.nanmin(grid1), np.nanmax(grid1), .01)
im = ax.contourf(y, z, grid1.T, cmap="nipy_spectral", levels = levels, extend="both")
ax.text(-60.5, 1900, "(8.%d) {}º".format(lon) % (i), va="center", ha="left", fontsize=12, color='black')
#ax.set_title("{} at {} W at {} resolution".format(title,lon, 1 / resolution))
#ax.set_title("at {}".format(lon))
cbticks = np.arange(dmin,dmax,((dmax-dmin)/9))
cb = fig.colorbar(im, pad=0.025, format='%.1f', extend='both',ticks=cbticks)
cb_ticklabel_size = 14 # Adjust as appropriate.
cb.ax.tick_params(labelsize=cb_ticklabel_size)
label_size = 16 # Adjust as appropriate.
#ax.tick_params(axis="both", direction='inout', which='major')
if i<3:
ax.tick_params(labelbottom=False)
ax.tick_params(labeltop=False, direction='inout')
else:
ax.tick_params(labelbottom=True)
ax.set_xlabel("° Latitude", fontsize=label_size)
if (i % 2) == 0:
ax.tick_params(labelleft=False)
cb.set_label(cbtitle, rotation=270, fontsize="16", labelpad=12.5)
else:
ax.set_ylabel("Depth (dbar)", fontsize=label_size)
contour_levels = np.linspace(tmin, tmax, num=10) # Adjust as appropriate.
print(contour_levels)
CS = ax.contour(y, z, grid2.T, #10,
colors='k', vmin=tmin, vmax=tmax, levels=contour_levels) # not sure if negative contours will be dashed by default
ax.clabel(CS, fontsize=9, inline=1, fmt='%1.2f')
#cb = fig.colorbar(im, CS, pad=0.025, format='%.2f', extend='both',ticks=cbticks)
db = oceansdb.ETOPO()
y = np.arange(-60, -75, -.01)
h = -db['topography'].extract(lat=y, lon=long)["height"]
#ax.plot(y, h)
ax.fill_between(y, 5000, h, color='black')
ax.set_ylim(2000, 0)
ax.set_xlim(-60, -74.5)
fig = plt.figure(figsize=(25,10))
fig.subplots_adjust(hspace=0.05, wspace=0.005)
for i in range(1, 5):
a = z_grid[i]
long = longit[i]
print(a,long)
xsection_dens = dens[:,a,:]
xsection_dens = interp_nans(xsection_dens)
xsection_temp = combined_temp_grid_withdepth[:,a,:]
xsection_temp = interp_nans(xsection_temp)
xsection_sal = argo_sal_grid_withdepth[:,a,:]
xsection_sal = interp_nans(xsection_sal)
b,c,d,e,f,g = [np.nanmin(xsection_temp),np.nanmax(xsection_temp),
np.nanmin(xsection_sal),np.nanmax(xsection_sal),
np.nanmin(xsection_dens),np.nanmax(xsection_dens)]
print(b,c,d,e,f,g)
plot_cross_slope_dens(xsection_dens, xsection_temp,
lon = long, i = i, fig = fig,
tmin=b, tmax=c, smin=d, smax=e, dmin=f, dmax=g)
#fig.suptitle("Cross section of {} at {} resolution".format(title, 1 / resolution), fontsize="20")
plt.show()
0 180 -1.7390999794006348 3.3632362650500403 33.731198120117185 34.73633130391439 1026.9056454753093 1027.8418129262836 [-1.73909998 -1.17217373 -0.60524748 -0.03832123 0.52860502 1.09553127 1.66245752 2.22938377 2.79631002 3.36323627] 90 -90 -1.8242663070559502 3.86983315149943 33.45185068491343 34.73500061035156 1026.667654198818 1027.8249655817285 [-1.82426631 -1.19158859 -0.55891087 0.07376685 0.70644456 1.33912228 1.9718 2.60447772 3.23715543 3.86983315] 180 0.0 -1.7166638456541916 1.002409734641369 33.870741965553975 34.69793701171875 1027.2289829270503 1027.8443220775991 [-1.71666385 -1.41454456 -1.11242527 -0.81030599 -0.5081867 -0.20606741 0.09605187 0.39817116 0.70029045 1.00240973] 270 90 -1.8788625001907349 1.8338245834623064 33.369389257123395 34.735844135284424 1027.0517573601667 1027.8411490738727 [-1.8788625 -1.46634171 -1.05382093 -0.64130014 -0.22877935 0.18374144 0.59626222 1.00878301 1.4213038 1.83382458]
parameter = title = "Density (kg/m3)"
plot_cross_slope(dens[:,2,:], title=title, cbtitle=parameter, resolution=resolution)
#plot_cross_slope(interp_nans(dens[:,2,:]), title="Interpolated density at 179° W", cbtitle="Density", resolution=resolution)
#Cell added on 21/10/2019 - still testing
# adding function for TS Diagram
dens = dens0(argo_sal_grid_withdepth, argo_temp_grid_withdepth)
# print(combined_temp_grid_withdepth[0])
#shape of dens - (30,720,200) = (latitude, longitude, vertical layers)
cmap = plt.get_cmap('nipy_spectral');cmap.set_bad(color='white')
fig = plt.figure(4,figsize=(15,15)) #TS-diagram of Argo profiles in the Ross Sea Slice
#plt.scatter(argo_sal_grid_withdepth[0:10],argo_temp_grid_withdepth[0:10],s = 2,c=dens[0:10],cmap = cmap, edgecolors='none',alpha=0.8)
plt.scatter(argo_sal_grid_withdepth[5:14,:,0:9],argo_temp_grid_withdepth[5:14,:,0:9],s = 2,c=dens[5:14,:,0:9],cmap = cmap, edgecolors='none',alpha=0.8)
plt.xlim(np.nanmin(argo_temp_grid_withdepth[5:14,:,0:9]), np.nanmax(argo_temp_grid_withdepth[5:14,:,0:9]))
plt.xlim(np.nanmin(argo_sal_grid_withdepth[5:14,:,0:9]), np.nanmax(argo_sal_grid_withdepth[5:14,:,0:9]))
#plt.xlim(33.8,35);plt.ylim(-2.5,5)
plt.clb = plt.colorbar();plt.clb.ax.set_title('Density (kg/m$^3$)')
plt.xlabel('Pratical salinity',fontsize = 18);
plt.ylabel('Temperature ($^oC$)',fontsize = 18)
#ax.set_xlabel('Pratical salinity',fontsize = 18);
#ax.set_ylabel('Temperature ($^oC$)',fontsize = 18)
plt.title('T-S diagram, Argo profiles',fontsize = 20)
# plt.axis('tight')
plt.show()
# for i in range(0, len(FILES)):
# read = Dataset(source + count_60[j][1] + '/profiles/' + FILES[i], mode = 'r')
# Lat = read.variables['LATITUDE'][:]
# Lon = read.variables['LONGITUDE'][:]
# if Lon[0] < -140, Lon[0] > 160:
# Temp = read.variables['TEMP'][:]
# Psal = read.variables['PSAL'][:]
# sigma = dens0(Psal[0],Temp[0])
# scatter(Psal[0],Temp[0],s = 2,c=sigma,cmap = cmap, edgecolors='none',alpha=0.8)
# xlim(33.8,35);ylim(-2.5,5)
# clb = colorbar()#;clb.ax.set_title('Density (kg/m$^3$)')
# axis('tight');grid('on',alpha=0.2)
# # savepath = os.path.join('C:/Users/fvie285/Desktop/PhD_Project/Data/figures/Argo_South_of_50/aoml_figures/TS_diagram_RSS.png')
# # savefig(savepath, bbox_inches='tight', dpi=300)
# show()
fig, ax = plt.subplots(1, 1, figsize=(20,10))
print(argo_sal_grid_withdepth.shape)
seal_counts = np.sum(~np.isnan(seal_sal_grid_withdepth), axis=2)
argo_counts = np.sum(~np.isnan(argo_sal_grid_withdepth), axis=2)
print(np.argwhere(seal_counts + argo_counts >= 300))
depth = np.arange(0, 2000, 10)
ax.set_xlabel("Salinity")
ax.set_ylabel("Depth (dbar)")
ax.set_title("Salinity profile")
ax.set_ylim(2000, 0)
from scipy.signal import savgol_filter
def fill_nan(A):
'''
interpolate to fill nan values
'''
inds = np.arange(A.shape[0])
good = np.where(np.isfinite(A))
f = interpolate.interp1d(inds[good], A[good],bounds_error=False)
B = np.where(np.isfinite(A),A,f(inds))
return B
ax.plot(argo_sal_grid_withdepth[0, 0, :], depth)
ax.plot(seal_sal_grid_withdepth[0, 0, :], depth)
#ax.plot(savgol_filter(fill_nan(argo_sal_grid_withdepth[11, 466, :]), 51, 3), depth)
#ax.plot(savgol_filter(fill_nan(seal_sal_grid_withdepth[11, 466, :]), 51, 3), depth)
ax.legend(["argo", "seal"])
fig, ax = plt.subplots(1, 1, figsize=(10,20))
print(argo_sal_grid_withdepth.shape)
depth = np.arange(0, 2000, 10)
ax.set_xlabel("Salinity")
ax.set_ylabel("Depth (dbar)")
ax.set_ylim(2000, 0)
long = 2
long_label = (long/2)-180
lat_ini = "60ºS"
lat_end = "74.5ºS"
ax.set_title(a_title + "\nSalinity mean profile at {}º, from {} to {} latitude".format(long_label,lat_ini,lat_end))
for dim_lat in range(0,29):
try:
filtered = savgol_filter(fill_nan(argo_sal_grid_withdepth[dim_lat, long, :]), 51, 3)
ax.plot(filtered, depth)
except:
pass
#ax.legend(["argo[dim_lat]"])
#ax.plot(seal_sal_grid_withdepth[4, 715, :], depth)
#ax.legend(["argo", "seal"])
#print(argo_lats)
#print(argo_lats.shape)
#print(argo_lons.shape)
argo_temp_grid_withdepth_summer = np.load("argo_temp_grid_withdepth_Summer.npy")
argo_temp_grid_withdepth_winter = np.load("argo_temp_grid_withdepth_Winter.npy")
seal_temp_grid_withdepth_summmer = np.load("seal_temp_grid_withdepth_Summer.npy")
seal_temp_grid_withdepth_winter = np.load("seal_temp_grid_withdepth_Winter.npy")
combined_temp_grid_withdepth_summer = np.load("combined_temp_grid_withdepth_Summer.npy")
combined_temp_grid_withdepth_winter = np.load("combined_temp_grid_withdepth_Winter.npy")
argo_sal_grid_withdepth_summer = np.load("argo_sal_grid_withdepth_Summer.npy")
argo_sal_grid_withdepth_winter = np.load("argo_sal_grid_withdepth_Winter.npy")
seal_sal_grid_withdepth_summer = np.load("seal_sal_grid_withdepth_Summer.npy")
seal_sal_grid_withdepth_winter = np.load("seal_sal_grid_withdepth_Winter.npy")
combined_sal_grid_withdepth_summer = np.load("combined_sal_grid_withdepth_Summer.npy")
combined_sal_grid_withdepth_winter = np.load("combined_sal_grid_withdepth_Winter.npy")
for season in ["Summer", "Winter"]:
argo_temp = np.load("argo_temp_grid_withdepth_{}.npy".format(season))
seal_temp = np.load("seal_temp_grid_withdepth_{}.npy".format(season))
fig, ax = plt.subplots(1, 1, figsize=(10,20))
depth = np.arange(0, 2000, 10)
ax.set_xlabel("{} temperature °C".format(season))
ax.set_ylabel("Depth (dbar)")
ax.set_title("{} temperature °C profile".format(season))
ax.set_ylim(2000, 0)
ax.plot(argo_temp[6, 233, :], depth)
ax.plot(seal_temp[6, 233, :], depth)
ax.legend(["argo", "seal"])
print(argo_temp[8,238])
resolution = 2
years = [2007, 2011, 2015]
try:
grids_per_year = np.load("grids_per_year.npy", allow_pickle=True).item()
except:
grids_per_year = {}
for y in years:
for s in ["Summer", "Winter"]:
start = datetime(y,1,1)
end = datetime(y + 4,1,1)
argo_temp_grid_withdepth_yearfiltered = load_netcdf_grid(with_depth=True, resolution = resolution, filter_start = start, filter_end = end, filter_season = s)
seal_temp_grid_withdepth_yearfiltered = load_netcdf_grid("seal", with_depth=True, resolution = resolution, filter_start = start, filter_end = end, filter_season = s)
combined_temp_grid_withdepth_yearfiltered = np.nanmean((argo_temp_grid_withdepth_yearfiltered, seal_temp_grid_withdepth_yearfiltered), axis=0)
if y not in grids_per_year:
grids_per_year[y] = {
"Winter": {},
"Summer": {}
}
grids_per_year[y][s] = {
"argo": argo_temp_grid_withdepth_yearfiltered,
"seal": seal_temp_grid_withdepth_yearfiltered,
"combined": combined_temp_grid_withdepth_yearfiltered
}
np.save("grids_per_year", grids_per_year)
plt.rcParams['figure.facecolor']='white'
db = 50
for i, y in enumerate(years):
for s in ["Summer", "Winter"]:
argo_temp_grid_withdepth_yearfiltered = grids_per_year[y][s]["argo"]
seal_temp_grid_withdepth_yearfiltered = grids_per_year[y][s]["seal"]
combined_temp_grid_withdepth_yearfiltered = grids_per_year[y][s]["combined"]
title = f"at {db * 10}dbar from the year {y}-{y + 4} in {s}"
plot_grid(argo_temp_grid_withdepth_yearfiltered[:,:,db], title=a_title + "\n" + title, resolution=resolution, vmin=tmin, vmax=tmax)
plot_grid(seal_temp_grid_withdepth_yearfiltered[:,:,db], title=s_title + "\n" + title, resolution=resolution, vmin=tmin, vmax=tmax)
plot_grid(combined_temp_grid_withdepth_yearfiltered[:,:,db], title=comb_title + "\n" + title, resolution=resolution, vmin=tmin, vmax=tmax)
interp = interp_nans(combined_temp_grid_withdepth_yearfiltered[:,:,db])
plot_grid(interp, title=interp_title + "\n" + title, resolution=resolution, vmin=tmin, vmax=tmax)
if i > 0:
delta = argo_temp_grid_withdepth_yearfiltered[:,:,db] - grids_per_year[years[i - 1]][s]["argo"][:,:,db]
print(np.nanmin(delta), np.nanmean(delta), np.nanmax(delta))
plot_grid(delta, title=comb_title + "\n" + f"Delta from {years[i - 1]} - {y} in {s} at {db}dbar", resolution=resolution)
print(grids_per_year.keys())
delta = grids_per_year[2015]["Winter"]["argo"][:,:,db] - grids_per_year[2011]["Winter"]["argo"][:,:,db]
print(np.nanmax(delta))
print(np.where(delta == np.nanmax(delta)))
delta[10,636]
print(grids_per_year[2015]["Winter"]["argo"][10,636,db])
print(grids_per_year[2011]["Winter"]["argo"][10,636,db])
print(10/2, 636 /2 - 180)
mask = ((argo["lat"] < -64.75) & (argo["lat"] > -65.25) & (argo["lon"] > 137.25) & (argo["lon"] < 138.75))
for i in np.flatnonzero(mask):
pres = argo["pressure"][i] / 10
pres = np.floor(argo["pressure"][i] / 10).astype(int)
depthmask = pres == 50
temps = np.copy(argo["temperature"][i][depthmask])
print(f'i={i},lat={argo["lat"][i]},lon={argo["lon"][i]},dt={argo["datetime"][i]},temps={np.around(temps,2)},mean={np.around(np.nanmean(temps), 2)}')
mask = ((argo["lat"] < -64.75) & (argo["lat"] > -65.25) & (argo["lon"] > 137.25) & (argo["lon"] < 138.75))
fig, ax = plt.subplots(1, 1, figsize=(10,20))
ax.set_xlabel("Salinity")
ax.set_ylabel("Depth (dbar)")
ax.set_ylim(2000, 0)
def smooth(y, box_pts):
box = np.ones(box_pts)/box_pts
y_smooth = np.convolve(y, box, mode='same')
return y_smooth
for i in np.flatnonzero(mask):
sal = argo["salinity"][i]
pressure = argo["pressure"][i]
#sal = savgol_filter(sal, 69, 1)
ax.plot(sal, pressure)